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A  parameter  optimization  algorithm  is  developed  based 
on  the  gradient  projection  concept.  The  algorithm  is  applied  to 
the  aeroelastic  optimization  of  a  semi-infinite  panel  in  super¬ 
sonic  flow  and  to  the  aeroelastic-stress  optimization  of  a  delta¬ 
wing  in  supersonic  flow.  The  total  weight  is  taken  as  the  ob¬ 
jective  function.  Constraints  are  placed  on  flutter  speed  and 
minimum  thickness  for  the  panel  problem  and  on  flutter  speed, 
stress  and  minimum  thickness  for  the  delta  wing  problem. 

Sandwich  plate  structural  idealization  is  used  for  both 
the  panel  and  the  delta  wing  problem.  During  the  optimization 
the  core  thickness  is  held  constant  while  the  cover  skin  thick¬ 
ness  is  varied.  Finite  element  representation  is  used  to  model 
the  structures.  Constant  thickness  and  tapered  sandwich  beam 
elements  are  used  for  the  panel  problem.  A  high  precision 
triangular  plate  bending  element,  developed  by  Cowper,  Kosko, 
Lindberg  and  Olson,  is  used  for  the  delta-wing  problem.  The 
element  stiffness  and  mass  matrices  are  modified  for  linear 
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cover  skin  thickness  variation.  The  element  or  the  nodal  point 
cover  skin  thicknesses  are  used  as  design  variables  (optimi¬ 
zation  parameters) . 

Two-dimensional,  quasi-steady  aerodynamic  theory  is 
used  to  obtain  the  aerodynamic  forces .  The  aerodynamic  damping 
is  included.  The  same  finite  element  functions  are  used  in 
deriving  the  finite  element  mass,  stiffness,  aerodynamic  and 
damping  matrices,  i.e.  they  are  consistent. 

In  order  to  obtain  the  flutter  constraint,  system 
equations  of  motion  are  derived.  These  equations  form  a  qua¬ 
dratic  eigenvalue  problem  which  can  be  reduced  to  a  linear 
eigenvalue  problem  where  the  positive  real  parts  of  the  eigen¬ 
values  indicate  flutter.  The  optimization  procedure  required 
calculation  of  constraint  gradients.  Using  analytical  expres¬ 
sions,  derivatives  of  the  eigenvalues  with  respect  to  design 
variables  are  calculated  efficiently. 

For  stress  constraints  only  static  loading  is  con¬ 
sidered.  Von  Mises  yield  criterion  is  used  to  derive  the 
equations  for  stress  constraints .  Stress  constraints  are 
placed  at  nodal  points. 

Due  to  nonlinearity  of  the  constraints  the  optimi¬ 
zation  procedure  incorporates  constraint  tolerances  and  pro¬ 
vision  for  returning  to  the  constraints  when  they  are  violated. 

Different  finite  element  combinations  are  studied 
for  the  panel  problem.  Using  six  tapered  elements  across  the 
span,  the  effect  of  damping  on  the  optimum  shape  is  investi¬ 
gated  . 
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For  the  delta  wing  problem  a  3  x  3  grid  of  triangu¬ 
lar  bending  elements  is  used.  The  natural  boundary  conditions 
are  placed  on  curvatures  for  the  free  edges  as  well  as  geo¬ 
metric  boundary  conditions  for  the  clamped  edge.  The  effect 
of  constraints  on  the  optimum  shape  is  studied  by  optimizing 
the  structure  first  using  only  flutter  and  thickness  con¬ 
straints,  then  using  only  stress  and  thickness  constraints, 
and  finally  using  all  three  sets  of  constraints.  Total  weight 
reductions  of  47%  to  63%  are  obtained  depending  on  the  combi¬ 
nations  of  constraints. 
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CHAPTER  1 


INTRODUCTION 

1.1  Introduction  to  the  Subject  of  Aeroelastic  Optimization 

The  subject  of  this  thesis  is  aeroelastic  optimiza¬ 
tion,  which  falls  under  the  general  subject  of  structural 
optimization  with  eigenvalue  (dynamic)  constraints. 

Although  there  are  few  earlier  papers  (mainly  ref¬ 
erence  25) ,  the  keen  interest  in  this  subject  started  in  the 
late  1960's.  Papers  published  on  aeroelastic  optimization 
since  then  present  different  approaches  to  the  different  aspects 
of  the  problem,  namely  the  structural  formulation,  the  aero¬ 
dynamic  formulation,  the  strength  formulation,  the  formulation 
of  the  optimization  problem  and  the  solution  of  the  optimization 
problem. 

Generally  the  approach  to  the  structural  formulation 
falls  into  two  groups; 

1)  Continuous  formulation 

2)  Discrete  formulation 

Several  theories  are  available  for  the  aerodynamic 
formulation; 

1)  Strip  theory 

2)  Piston  theory 

3)  Quasi-steady  aerodynamic  theory 

4)  Three  dimensional  aerodynamic  theory 

Formulation  of  the  optimization  problem  follows  from 
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the  structural  formulation; 

1)  Control  theory  formulation  is  applied  when  the  con¬ 
tinuous  structural  formulation  is  employed. 

2)  Parameter  optimization  formulation  follows  the  dis¬ 
crete  structural  formulation. 

There  exist  several  approaches  to  the  solution  of  the 
control  theory  formulation.  Among  them  are  the  following; 

1)  Transition  Matrix  Method 

2)  Gradient  projection  method 

The  solution  methods  for  the  parameter  optimization 
problem  have  been  independently  developed  by  mathematicians  and 
numerical  analysts.  Reference  13  is  a  good  source  for  the 
numerical  techniques  available  for  structural  optimization. 

There  are  two  different  approaches  depending  on  the  way  that 
the  eigenvalue  constraint  is  handled; 

1)  If  the  eigenvalue  constraint  is  assumed  to  be  an 
equality  constraint,  then  the  problem  is  reduced  to  a  set  of 
non-linear  algebraic  equations  by  introducing  Lagrange  multi¬ 
pliers.  Several  iterative  numerical  techniques  exist  for  the 
solution  of  non-linear  algebraic  equations  such  as  the  Newton- 
Raphson  method  and  its  variations . 

2)  If  the  eigenvalue  constraint  is  assumed  to  be  an 
inequality  constraint,  then  the  problem  can  be  stated  as  a  non¬ 
linear  mathematical-programming  problem.  One  of  the  iterative 
techniques  based  on  gradient  information  can  be  used  for  the 
solution.  The  most  commonly  used  algorithms  in  the  aeroelastic 
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optimization  field  are; 

1)  Interior  penalty  function  method 

2)  Usable-feasible  directions  method  of  Zoutendijk. 

3)  Gradient  projection  method 

The  predominant  constraint  in  aeroelastic  optimization 
is  flutter.  Few  of  the  works  published  on  aeroelastic  optimi¬ 
zation  also  consider  other  constraint  such  as  deflections  and 
stress.  The  process  of  handling  such  constraints  has  already 
been  well  established  in  the  field  of  structural  optimization. 

In  the  following  some  of  the  important  works  on  the 
aeroelastic  optimization  are  briefly  discussed: 

In  reference  4  Turner  has  extended  his  earlier  work 
on  structural  optimization  with  constraints  on  natural  frequen¬ 
cies  (Ref.  26)  to  the  optimization  of  structures  with  flutter 
constraints.  He  uses  the  discrete  approach  but  handles  the 
flutter  constraint  as  an  equality  constraint.  He  obtains  a  set 
of  non-linear  algebraic  equations  by  using  complex  Lagrange 
multipliers.  The  equations  are  solved  using  the  Newton- Raphson 
method.  Two  examples  are  presented:  1)  Panel  Optimization 
2)  Cantilever  Wing  Optimization.  His  method  can  be  used  with 
any  aerodynamic  theory.  He  uses  quasi-steady  aerodynamic  theory 
for  the  panel  problem  and  strip  theory  for  the  cantilevered 
wing.  The  number  of  design  variables  is  kept  rather  small 
(3  elements  with  constant  mass  and  stiffness  properties  are 
used  for  both  problems) ,  presumably  because  of  computational 
effort  involved  in  solving  non-linear  equations.  No  other  con- 
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straints  are  considered. 

A  paper  by  Ashley  and  McIntosh  (Ref.  1)  presents  re¬ 
sults  using  the  continuous  formulation.  Several  examples  with 
eigenvalue  constraints  are  discussed,  among  which  are  the  opti¬ 
mization  of  a  two-dimensional  plate  and  of  a  cantilevered  lift¬ 
ing  surface  with  flutter  constraints.  Quasi-steady  aerodynamic 
theory  with  zero  aerodynamic  damping  is  used.  Equations  are 
derived  using  control  theory  formulation  but  left  unsolved  for 
the  two-dimensional  plate  problem.  An  analytical  solution  is 
given  by  inspection  for  the  cantilevered  wing.  Placing  con¬ 
straints  on  minimum  thickness  is  also  considered. 

Weisshaar,  in  reference  7,  made  significant  contri¬ 
butions  to  both  continuous  and  discrete  approaches.  Using  con¬ 
trol  theory  formulation,  he  obtains  simultaneous  non-linear 
differential  equations  and  solves  these  equations  using  the 
"transition  matrix"  procedure  (Ref.  27) .  Among  examples  with 
dynamic  constraints  he  presents  solutions  for  the  first  time  to 
the  optimum  skin  thickness  distribution  of  a  simply  supported 
sandwich  panel  in  supersonic  flow  with  different  thickness  con¬ 
straints  and  different  mass  ratios.  Quasi-steady  aerodynamic 
theory  is  employed.  The  aerodynamic  damping  is  neglected.  He 
also  presents  some  results  using  finite  elements  and  a  first- 
order  gradient  technique  (based  on  Rubin's  work  in  reference  28) 
on  the  same  problem.  Flutter  constraint  gradients  are  calcu¬ 
lated  numerically. 

A  paper  by  Bhatia  and  Rudisill  (Ref.  9)  present  their 
work  on  the  optimization  of  complex  structures  with  flutter 
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constraints  using  finite  elements.  The  optimization  procedure 
consists  of  three  different  search  procedures,  the  first  two 
are  based  on  steepest-descent  principle;  one  to  increase  the 
flutter  speed,  the  other  to  decrease  the  structural  weight. 

The  third  search  procedure  is  based  on  gradient  projection  con¬ 
cept;  it  holds  the  structural  weight  constant  while  maximizing 
flutter  speed.  The  main  contribution  is  the  analytical  expres¬ 
sions  for  the  calculation  of  flutter  constraint  gradients.  A 
cantilevered  box  beam  is  optimized  using  three  dimensional 
aerodynamic  theory.  With  only  flutter  and  thickness  constraints 
approximately  40%  total  weight  reduction  is  obtained  after  22 
design  cycles. 

A  paper  by  Craig  (Ref.  8)  present  results  for  flutter 
optimization  of  a  simply  supported  sandwich  panel  using  finite 
elements.  Consistent  element  aerodynamic  matrices  are  derived 
using  quasi-steady  theory  with  zero  damping.  Two  different 
optimization  modes  are  employed  both  based  on  the  gradient  pro¬ 
jection  concept.  The  weight  minimization  mode  minimizes  the 
weight  while  holding  the  flutter  parameter  constant.  The 
Lambda-modification  mode  increases  the  flutter  parameter  while 
holding  the  weight  constant.  Optimum  shapes  are  given  for  five 
and  nine  constant  thickness  beam  elements  across  the  span. 

Rao  in  reference  19  presents  a  procedure  for  optimum 
design  of  aircraft  wings  to  satisfy  strength,  stability,  fre¬ 
quency  and  flutter  constraints.  Constant  stress  membrane  ele¬ 
ments  and  rectangular  shear  panels  are  used  for  structural 
idealization.  Piston  theory  with  zero  damping  is  used  to  calcu- 
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late  aerodynamic  forces.  An  interior  penalty  function  method 
is  applied  to  the  optimization  problem.  Eigenvalue  gradients 
are  calculated  using  exact  expressions.  Optimum  shapes  are 
obtained  for  several  supersonic  aircraft  wings  using  4  consecu¬ 
tive  penalty  parameters  and  25  to  30  design  cycles. 

Miura,  in  reference  20,  discusses  a  procedure  to  de¬ 
termine  the  optimum  configuration  of  a  lifting  surface  with 
aeroelastic  constraints.  Sandwich  plate  structural  idealization 
is  used.  Ritz  type  displacement  (modal)  method  is  applied  to 
derive  system  matrices.  Aerodynamic  forces  are  obtained  using 
piston  theory  with  negligible  damping.  The  usable-feasible 
directions  method  of  Zoutendijk  is  used  to  solve  the  optimiza¬ 
tion  problem.  Constraints  are  also  placed  on  stresses  and  de¬ 
flections.  The  paper  by  Fox,  Miura  and  Rao  (Ref.  10)  summaries 
the  above  works  by  Rao  and  Miura. 

A  report  put  out  by  the  Langley  Research  Center  (Ref. 
29)  gives  the  optimum  thickness  distributions  of  a  delta  wing 
for  stress  and  thickness  constraints,  for  flutter  and  thickness 
constraints  and  for  stress,  flutter  and  thickness  constraints. 
Sandwich  plate  idealization  and  piston  theory  are  used.  A 
finite-difference  method  is  employed  for  discretization.  An 
interior  penalty  function  approach  is  chosen  for  the  solution 
of  the  optimization  problem. 

Gwin,  in  reference  21,  presents  computer  routines  for 
the  minimum-weight  synthesis  of  a  lifting  surface  with  flutter 
constraint.  Three-dimensional  aerodynamic  theory  is  used. 
Different  finite  elements  can  be  used  to  represent  the  struc- 
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ture.  A  modified  version  of  Zoutendijk's  method  is  coded  to 
solve  the  parameter  optimization  problem.  System  matrices  are 
reduced  using  the  first  few  natural  vibration  modes  as  gen¬ 
eralized  coordinates.  Analytical  expressions  are  used  for 
flutter  constraint  gradients.  Several  subsonic  examples  are 
presented. 

Recently,  in  reference  30,  Pierson  presented  new  con¬ 
tributions  to  the  continuous  approach.  He  applied  what  he 
calls  a  gradient  projection  optimal  control  method  to  control 
theory  formulations.  Results  are  presented  for  one  and  two 
dimensional  solid  and  sandwich  panels  with  different  boundary 
conditions.  Quasi-steady  aerodynamic  theory  is  used.  A  dis¬ 
cussion  is  given  to  handle  aerodynamic  damping  but  the  pre¬ 
sented  results  are  for  zero  damping. 

1.2  Introduction  to  the  Present  Study: 

In  light  of  the  above  discussion  the  present  work  can 
be  briefly  described  as  follows: 

The  discrete  approach  is  adopted  for  the  structural 
formulation.  Sandwich  plate  and  finite  element  structural 
idealizations  are  used  to  represent  the  structure.  Thickness 
parameters  either  associated  with  finite  elements  or  with  nodal 
points  are  used  as  design  variables.  The  optimization  problem 
is  formulated  as  a  non-linear  mathematical-programming  problem 
where  the  objective  function  is  the  total  weight  of  the  struc¬ 
ture  which  is  linear  in  the  design  variables.  Flutter,  stress 
and  thickness  constraints  are  considered.  Quasi-steady  aero- 


dynamic  theory  with  damping  is  used  to  calculate  aerodynamic 
forces.  However  other  aerodynamic  theories  can  also  be  incor¬ 
porated  within  the  overall  design  procedure. 

The  consistent  mass,  stiffness,  aerodynamic  and  damp¬ 
ing  matrices  are  derived  for  elements  and  for  the  system.  The 
flutter  constraint  is  represented  as  a  quadratic  eigenvalue 
problem  in  terms  of  these  matrices.  This  is  reduced  to  a  linear 
eigenvalue  problem  where  the  eigenvalues  are  the  complex  fre¬ 
quencies  and  positive  real  parts  of  the  eigenvalues  indicate 
flutter.  It  is  possible  to  detect  flutter  phenomena  involving 
modes  other  than  the  two  initial  modes  during  the  optimization. 
This  way  of  handling  the  flutter  constraint  is  unique  to  this 
study.  In  previous  works  on  aeroelastic  optimization  the 
flutter  constraint  was  placed  on  the  flutter  speed.  That  makes 
it  very  hard  to  predict  flutter  involving  modes  other  than  the 
initial  two  modes.  Here  the  flutter  speed  is  determined  for  the 
initial  design  and  treated  as  a  constant  parameter  during  the 
optimization. 

The  stress  constraints  are  only  considered  for  uni¬ 
form  static  loading,  employed  as  a  first  approximation  to  the 
aeroelastic  loading.  This  is  consistent  with  the  quasi-steady 
aerodynamic  theory  approximation.  Von  Mises  yield  criterion 
is  used  to  obtain  stress  constraint  equations.  Constraints  are 
placed  at  finite  element  nodal  points.  Thickness  constraints 
are  directly  applied  on  design  variables. 

An  algorithm  based  on  the  gradient  projection  concept 
is  developed  to  solve  the  optimization  problem.  Due  to  non- 
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linearity  of  the  constraints,  the  algorithm  incorporates  con¬ 
straint  tolerances.  The  design  change  vector  is  composed  of 
two  orthogonal  vectors,  one  in  the  projected  gradient  direction 
to  reduce  the  objective  function;  the  other  in  the  orthogonal 
direction  for  returning  to  the  violated  constraints. 

The  optimization  procedure  requires  calculation  of 
constraint  gradients.  Analytical  expressions  are  derived  to 
perform  these  calculations  efficiently. 

Two  example  problems  are  considered:  A  semi-infinite 
panel  in  supersonic  flow  and  a  cantilevered  delta  wing  in  super¬ 
sonic  flow.  Different  finite  element  representations  are  tried 
for  the  panel  problem.  Results  are  compared  with  Weisshaar's 
(Ref.  7)  and  Craig's  (Ref.  9)  results.  Six  tapered  sandwich 
beam  elements  are  used  across  the  panel  span  to  study  the 
effect  of  damping  on  the  optimum  shape. 

For  the  delta  wing  problem  high  precision  triangular 
plate  bending  elements  (Ref.  22)  are  used.  Natural  frequencies 
and  flutter  boundaries  are  checked  against  results  given  by 
Olson  and  others  (Refs.  6  and  22) .  The  minimum-weight  designs 
are  obtained  for  flutter  and  thickness  constraints,  for  stress 
and  thickness  constraints  and  for  flutter,  stress  and  thickness 
constraints. 

The  main  objectives  of  this  study  can  be  summarized 
as  follows : 

1)  Development  of  a  parameter  optimization  algorithm 
for  the  minimum-weight  design  of  lifting  surfaces;  an  algorithm 
capable  of  giving  an  approximate  design  using  a  small  number 
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of  iterations. 

2)  Development  of  equations  for  handling  the  flutter 
constraint  such  that  is  is  possible  to  predict  flutter  involv¬ 
ing  higher  modes  during  the  optimization. 

3)  Development  of  analytical  expressions  for  the  con¬ 
straint  gradients  so  that  these  can  be  calculated  with  minimum 
time . 

4)  Study  the  effect  of  aerodynamic  and  structural 
damping  on  the  minumum-weight  design. 

5)  Study  the  effect  of  different  combinations  of 
flutter,  stress  and  thickness  constraints  on  the  minimum-weight 
design. 

1.3  Organization  of  the  Thesis; 

In  Chapter  2  the  optimization  formulation  is  pre¬ 
sented.  Several  popular  algorithms  for  the  solution  of  the  non¬ 
linear  mathematical-programming  problem  are  briefly  discussed. 

Chapter  3  is  for  structural  formulation:  Sandwich 
idealization  is  presented.  Equations  for  stresses  and  bending 
moments  in  terms  of  transverse  displacement  second  derivatives 
are  given.  Finite  element  representation  is  briefly  explained. 

Equations  for  aerodynamic  forces  are  derived  in  Chap¬ 
ter  4.  The  derivation  of  element  aerodynamic  and  damping 
matrices  are  discussed. 

In  Chapter  5  equations  expressing  flutter,  stress  and 
thickness  constraints  in  terms  of  design  variables  are  derived. 
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Chapter  6  gives  the  derivation  of  equations  for  the 
gradient  projection  algorithm.  Figures  explaining  different 
design  steps  are  also  presented. 

The  analytical  expressions  for  flutter,  stress  and 
thickness  constraint  gradients  are  developed  in  Chapter  7. 

Chapter  8  and  Chapter  9  present  the  results  for  the 
semi-infinite  panel  optimization  and  the  cantilevered  delta 
wing  optimization  respectively. 

In  Chapter  10  several  conclusions  and  suggestions  for 
further  research  are  listed. 

Appendices  A,  B  and  C  are  for  the  derivation  of  finite 
element  matrices  for  the  different  finite  elements  used  in  the 


examples . 


CHAPTER  2 


FORMULATION  OF  THE  OPTIMIZATION  PROBLEM 


Design  of  aircraft  structures  can  be  cast  as  a  non¬ 
linear  mathematical-programming  problem  of  the  form: 

Find  p 

such  that  F(p)  -»•  minimum 

subject  to  Cj(p)  £  0;  j  =  1,2,...,  m 

Here  p  ,  the  vector  of  optimization  parameters, 
should  contain  well  chosen  design  variables  such  that  when  p 
is  determined,  the  structure  is  essentially  determined.  F(p), 
the  merit  (objective)  function,  should  be  so  chosen  that  its 
minimization  would  drive  us  to  a  most  desirable  structural 
configuration,  mainly  a  most  economic  one.  (p) , 
j  =  1,  2,  ...,  m,  the  constraints,  should  include  all  the  cri¬ 
teria  the  structure  has  to  satisfy  for  its  usability. 

Most  established  algorithms  for  the  solution  of  the 
non-linear  programming  problem  require  the  calculations  of  the 
merit  function  gradient; 


VF  (n } T  -  f8F  8F  8F  \ 

VF  ~  “  {8px  '  8p2  '**•'  3p^}(P) 


(2.1) 


and  the  constraint  gradients; 


VC.  (p) 


sp-l  '  9p2 


9C  . 

*'  9^}(p)  '  3  =  1,2,. ..,m 

11  ~ 


(2.2) 
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since  along  these  directions  the  changes  in  F(p)  and  Cj (p)  , 

j  =  1,2,...,  m,  are  maximum.  For  this  reason  the  functional 

relationship  between  F(p),  C.(p),  j  =  1,2,...,  m,  and  p  should 

3  ~ 

be  mathematically  known  to  us. 

The  following  definitions  are  in  order  for  the  clarity 
of  the  rest  of  this  chapter: 

1)  Any  design  p^  is  said  to  be  feasible  if  it  satis¬ 
fies  the  relationships; 

C. (p  )  £  0;  j  =  1,2,...,  m  (2.3) 

J  ^4 

2)  Any  design  p  is  said  to  be  usable  with  respect  to 

~q 

a  previous  design  p  if  the  relationship; 

~  M.  ■** 

F  (Pg)  £  F(Pq_i)  (2  • 4) 

is  valid. 

3)  A  usable-feasible  design  change  satisfies  both  of 
these  conditions  simultaneously. 

4)  A  constraint  is  said  to  be  "convex"  if,  any  move 
made  from  the  constraint  along  a  direction  orthogonal  to  the 
gradient  vector,  leads  to  a  design  in  the  infeasible  region. 

In  figure  2.1  some  of  the  concepts  described  above 
are  illustrated  on  a  2-dimensional  design  space  (n  =  2) . 

2.1  Design  Variables  (Optimization  Parameters) 

Generally,  sizing  parameters  of  the  structural  mem¬ 
bers  are  used  as  design  variables  in  structural  optimization. 

The  finite  element  representation  makes  the  choice  particularly 
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easy.  In  the  problems  considered  here  either  the  finite  element 
thicknesses  or  the  nodal  point  thicknesses  are  used  as  design 
variables . 

2.2  Merit  (Objective)  Function 

The  structural  weight  is  chosen  to  be  the  merit  func¬ 
tion  in  all  the  cases  considered  here.  This  choice  is  rather 
common  for  aeroelastic  optimization  since  both  the  initial  and 
operational  costs  of  aircraft  structures  are  closely  related  to 
the  structural  weight. 

As  shall  be  seen  in  the  example  problems  it  is  very 
easy  to  express  the  structural  weight  in  terms  of  either  the 
finite  element  thicknesses  or  the  nodal  point  thicknesses.  In 
all  cases  the  structural  weight  will  be  a  linear  function  of 
the  design  variables.  This  fact  makes  the  calculations  of 
F(p),  the  merit  function,  and  of  VF(p),  the  merit  function 
gradients  particularly  simple  and  time-saving. 

2.3  Constraints 

Most  aircraft  structures  are  required  to  satisfy  the 
following  types  of  design  criteria: 

1.  Combined  stresses  must  be  below  yield  for  all 
structural  members  and  for  all  specified  loading  conditions. 

2.  Deflections  under  all  specified  loading  condi¬ 
tions  must  be  within  allowable  limits. 

3.  Member  sizes  must  be  within  practical  limits. 

4.  Natural  frequencies  must  fall  within  allowable 
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bands . 

5.  The  structure  must  be  aerodynamically  stable, 

i.e.,  it  must  be  free  from  flutter  and  divergence. 

The  following  constraints,  taken  from  the  above  list, 
are  considered  in  the  applications  discussed  here: 

1.  Flutter:  The  flutter  speed  of  the  final  minimum- 
weight  design  must  be  greater  than  or  equal  to  the  flutter 
speed  of  the  initial  design. 

2.  Stress:  The  final  minimum-weight  design  must 
satisfy  the  same  yield  criteria  satisfied  by  the  initial  design 

3.  Thickness:  The  member  sizes  of  the  final  minimum 
weight  design  must  be  greater  than  or  equal  to  a  certain  per¬ 
centage  of  the  member  sizes  of  the  initial  design. 

In  the  examples  discussed  in  Chapter  8  and  Chapter  9, 
the  effect  of  different  combinations  of  constraints  on  the 
optimization  is  studied.  Flutter  +  thickness,  stress  +  thick¬ 
ness,  and  flutter  +  stress  +  thickness  combinations  are  con¬ 
sidered.  Chapter  5  contains  the  explicit  expressions  for  the 
constraints  as  functions  of  p.  Stress  and  flutter  constraints 
are  non-linear;  thickness  constraints  are  linear  constraints. 

2.4  Optimization  Algorithms 

Three  different  algorithms,  popular  in  structural 
optimization,  are  considered  to  drive  the  optimization  process 
from  pQ,  the  initial  design  variables,  to  p^,  the  final  (mini¬ 
mum-weight)  design  variables.  These  algorithms  are  briefly 
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discussed  below: 

1.  Usable-Feasible  Directions  Method  of  Zoutendijk  (Ref.  13, 
Sec.  7.2): 

The  idea  is  to  generate  a  sequence  of  feasible  designs 
Pq/  p^,...,  Pg,  Pq+i  such  that  for  all  q  the  inequality 


F(eq+i>  i  F'eq> 


is  satisfied.  Here  Pg+^  is  defined  as: 


)  ,  =  p  +  a  S. 

:q+l  ~q  q  ~q 


(2.5) 


(2.6) 


so  that  the  qth  step  of  the  process  involves  the  calculation  of 

,  ,  * 

the  direction  vector  Sg  and  the  step  length  ag  . 

C.(p  )  is  a  critical  constraint  at  p  if  C. (p  )  =  0. 

3  4  J  M 

For  such  critical  constraints  the  feasibility  of  is  assured 

if  S  satisfies  the  inequality: 

~q 

S  TVC. (p  )  <  0  ;  j  e  J  (2.7) 

~q  3  ~q  ° 


where  J  is  the  set  of  all  critical  constraints, 
c 

The  usability  of  S  ,  on  the  other  hand,  is  assured  if 

~q 

Sg  satisfies  the  inequality: 

S  TVF(p  )  <  0  (2.8) 

~q 

The  determination  of  S  satisfying  the  above  inequali- 

— 

ties  can  be  cast  as  an  auxiliary  optimization  problem  as  follows: 
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Maximize  8 

Subject  to: 

(1) 

?qTvcj(Pq>  +  V  1  0  ;  j  e  Jc 

(2.9) 

(2) 

s  tvf(p  )  +  8  <  0 

~q  ~q  — 

(2.10) 

(3) 

is  normalized 

~q 

Here 

0 .  is  a  pre-selected  non-negative 

parameter  to 

account  for  the  curvature  of  a  non-linear  convex  constraint. 

0j  =  0  is  used  for  a  linear  constraint. 

This  auxiliary  optimization  problem  is  linear,  pro¬ 
vided  that  a  linear  normalization  is  used  for  S  ,  such  as, 

~q 

-1  £  S  .  £  1.  Therefore  the  elements  of  S  and  6  can  be  easily 
calculated  by  using  a  linear  programming  algorithm  (simplex 
method).  If  the  auxiliary  problem  leads  to  8  >  0,  then  F(p) 
can  be  improved  within  the  feasible  region.  If,  however,  we 
obtain  8=0  then  it  can  be  shown  that  p  is  the  optimal  solu- 
tion. 

Once  the  direction  S  has  been  determined  as  above, 

determination  of  a*  for  non-linear  merit  functions  can  be  cast 

4 

as  a  one-dimensional  minimization  problem  (Ref.  13,  sec.  6.2.2); 

Find  a  such  that  F(p^  +  a  S )  =  F (a  )  is  minimum  provided 

q  q~q  q 

P  , i  =  p  +  a  S  is  in  the  feasible  region.  In  practice 
~qT_L  ~q  q~q 

F (a  )  is  usually  approximated  by  a  quadratic  function  and  deter- 
~1 

* 

mination  of  a  requires  calculation  of  F(a  )  for  two  different 

q  q 

values  of  a  other  than  a  =0.  If  the  gradients  of  the  merit 
4  4 

function  are  easily  obtained,  then  a  cubic  approximation  can  be 
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used  in  which  case  calculation  of  F(a  )  and  VF (a  )  for  one  value 
of  a  other  than  a  =  0  is  needed. 

q  q 

For  linear  merit  functions  the  only  consideration  for 
the  selection  of  a*  is  the  feasilibility  of  Pg+^*  Therefore 
a*  can  be  supplied  by  the  user  (Ref.  21). 

When  a  is  so  obtained,  the  constraints  are  checked. 

If  there  is  a  constraint  violation,  an  interpolation  procedure 
must  be  used  to  reduce  the  a  so  that  we  are  on  or  near  the 

q 

violated  constraint,  within  the  feasible  region  (Ref.  21) . 

One  way  to  increase  the  efficiency  of  the  algorithm 

is  to  include  the  near-critical  constraints,  -e  <  C. (p  )  <  0 

—  j  ~q  — 

into  the  set  of  critical  constraints.  Here  e  is  a  positive 
parameter  which  can  be  made  successively  smaller  when  small 
values  of  3  indicate  that  p  approaches  the  optimal  solution. 

It  should  be  noted  that  the  choice  of  scaling  coef¬ 


ficients  0j,  j  =  1,2,...,  m,  for  non-linear  constraints  is 
intuitive,  and  it  will  affect  the  convergence  speed  of  the 
algorithm  considerably.  Large  values  for  0j's  will  push  the 
iteration  away  from  the  constraints  and  most  likely  away  from 
the  optimum  point  too.  Small  values  for  0j's  on  the  other 
hand  will  be  ineffective  in  preventing  the  same  constraints 
from  being  violated  in  the  following  steps. 

The  usable-feasible  directions  algorithm  is  terminated 
when  either  3  or  the  change  in  F(p)  becomes  small. 
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2.  The  Gradient  Projection  Method  (Ref.  13,  Sec.  7.3): 

This  method  assumes  that  the  optimum  point  lies  on 
one  or  more  constraints.  Therefore  the  idea  is  to  follow  the 
constraints  as  closely  as  possible  while  reducing  the  merit 
function.  For  this  purpose  the  method  uses  the  projections  of 
the  steepest  descent  direction,  -VF(p  ),  into  the  manifold  de- 

''■'CJ 

fined  by  currently  active  constraints,  C. (p  )  =  0,  to  calcu- 

3  ~c[ 

late  the  vector  S  (same  meaning  as  in  equation  2.6).  Thus  S 

~q  ~q 

satisfies  the  relationships; 


and 


VF(p  )  < 


0 


Sq  VCi  V 


=  o  ;  j  e  J 


(2.11) 


(2.12) 


where  Jc  is  the  set  of  all  critical  constraints. 

It  is  possible  to  express  as: 

?q  “  -  ?q  VF<gq>  <2-13> 

where  is  the  projection  matrix.  P^  can  be  calculated  by 
using  the  projection  theorem  of  vector  calculus.  First  it  is 
convenient  to  introduce  the  rectangular  matrix  of  constraint 
gradients: 

(2.14) 


G 

~q 


VCl*£q*'  VC2^£q*,‘**'  VCm*£q* 


where  m  is  the  number  of  active  constraints,  so  that  we  can 
write  the  feasibility  condition  (equations  2.12)  concisely  as 
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G  S  =  0 

~q 


(2.15) 


The  projected  vector  P  VF(p  )  can  be  obtained  from 

~q  ~q 

VF(p  )  by  substracting  from  it  the  vector  G  v  ,  which  is  a 
linear  combination  of  the  constraint  gradients.  Now  we  can  re¬ 
write  the  feasibility  condition: 

SqT  [VF(-Pq,l  =  9qTpF(eq'  '  IqYql  =  9  (2'16) 


This  equation  leads  to: 


v  =  G  TG  G  T  VF(p  ) 

~q  ~q  ~q  ^q  ~q 


(2.17) 


T  T 

P  =  I  -  G  G  G  G 

~q  z  ^q  *q  *q  »q 


(2.18) 


We  should  note  that  the  matrix  [G  G  ]  is  an  m  x  m 

~  ~  M. 

square  matrix  and  can  generally  be  inverted  if  the  constraint 
gradients  are  linearly  independent. 

Once  S  is  so  determined  and  normalized  we  have  the 

~q 

usual  equation  for  qth  iteration  step; 


£q+l  ~q  +  aq  ~q 


(2.19) 


The  determination  of  a  can  again  be  cast  as  a  one- 
dimensional  minimization  problem  if  the  merit  function  is  non¬ 
linear,  if  the  merit  function  is  linear  it  can  be  given  a 
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priori.  For  either  case  a*  must  be  bounded  not  to  violate  any 
new  constraints.  Interpolation  can  be  used  to  reduce  a*  in 

q 


case  such  a  violation  occurs. 


As  formulated  above  the  gradient  projection  method  is 
very  powerful  for  linear  constraints  but  will  not  work  very 
well  for  non-linear  constraints.  When  we  are  on  a  convex  con¬ 
straint  any  move  performed  along  S  will  take  us  out  of  the 

~q 

feasible  region.  An  extra  step  has  to  be  performed  to  come 
back  to  the  feasible  region  (see  Chapter  6  for  a  detailed  dis¬ 
cussion)  .  For  such  cases,  the  choice  of  a*  becomes  very  impor- 

ic 

tant.  A  large  value  for  a  will  increase  the  difficulty  in 

Si 

coming  back  into  the  feasible  region.  On  the  other  hand,  a 

small  value  for  a  will  mean  that  too  many  iterations  will  be 

required  to  converge  to  the  optimum  solution. 

The  idea  of  including  the  near  critical  constraints, 

-e  <_  C.  (p  )  £  0  ,  into  the  set  of  critical  constraints,  to  in- 

crease  the  efficiency  of  the  algorithm  also  applies  for  the 

gradient  projection  method. 

The  convergence  of  the  algorithm  can  be  based  either 

on  the  amount  of  variation  of  the  objective  function  in  one 

iteration  step,  F(p  )  -F(p  . , ) ,  or  on  the  product  S  TVF(p  ) 

~q  ~q+i  ~q  ~q 

which  tends  to  get  smaller  as  the  optimum  point  is  approached. 


3.  The  Interior  Penalty  Function  Method  (Ref.  13,  Sec.  6.3): 

The  idea  is  first  to  convert  the  constrained  mini¬ 
mization  problem  into  an  unconstrained  minimization  problem  and 
then  solve  the  unconstrained  minimization  problem  by  one  of  the 
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well  developed  methods.  The  conversion  is  done  by  augmenting 
the  objective  function  with  a  penalty  term  which  is  small  at 
points  away  from  the  constraints  in  the  feasible  region.,  but 
which  gets  very  large  as  the  constraints  are  approached.  The 
most  commonly  used  form  is  the  following: 

m  . 

<M  P,r)  =  F  ( p )  -  r  l  x  ,  (2.20) 

j=l  <-j 

where  F(p)  is  to  be  minimized  over  all  p  satisfying 

C.(p)  <  0,  j  =  1,2,...,  m.  Here  the  penalty  parameter  r  is  a 
J  ~  — 

positive  number  which  must  be  made  successively  smaller  in  order 
to  come  as  close  to  the  constraints  as  possible. 

The  most  powerful  mehtod  so  far  developed  to  solve 
the  unconstrained  minimization  problem  is  the  Davidon-Fletcher- 
Powell  variable  metric  method  (Ref.  13,  sec.  6.2.6)  and  it  is 
usually  coupled  in  practice  with  the  interior  penalty  function 
approach  to  drive  the  optimization  for  decreasing  r  values. 

It  proceeds  as  follows: 


Bq+1 


=  p„  +  a 

tq  q  ~q 


(2.21) 


where  the  directional  vector  is  calculated  as: 

~q 

(2)  Sg  =  -  HqV<MPq,r)  (2.22) 

Here  -V<(>(pq,r)  is  the  well  known  steepest  descent 
direction.  The  matrix  H  is  a  modifier  which  takes  care  of 
the  phenomenon  called  zig  zagging  which  is  the  main  cause  of 
inefficiency  associated  with  steepest  descent  algorithms.  We 


can  write 
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H 

F  (H  ,  ,  S  .  , 
~q-l  ~q-l ' 

where  g  = 

~q 

V<f>  (pq,r) 

and  g  . 

~q-i 

=  V(Mpq_i,r) 

?q-l'  0tq-l)  (2.23) 

(2.24) 

(2.25) 


The  calculation  details  of  H  will  be  omitted  here 

~q 

but  may  be  found  in  reference  13.  It  can  be  shown  that  it  is  a 
positive  definite  approximation  to  the  matrix  of  second  partial 
derivatives  (Ref.  13) .  Ho  usually  is  taken  to  be  equal  to  the 
identity  matrix  I. 


(3)  Once  is  known,  is  calculated  by  solving 

the  one-dimensional  minimization  problem  (Ref.  13,  Sec.  6.2.2): 

Find  a  such  that  <f>(Pn  +  a  S  ,r)  e  <j)  (a  )  is  minimum,  provided 

p  ,  -i  =  p  +  a  S  is  in  the  feasible  region.  Since  <b(p,r) 
~q+l  ~q  q  ~q  T~ 

is  non-linear,  either  a  quadratic  or  a  cubic  approximation  can 
be  used  for  (j)  (ctg)  as  was  discussed  for  Zoutendijk's  method. 

Once  a*  is  determined  the  qth  step  is  completed. 

The  following  points  should  be  noted  about  the  in¬ 
terior  penalty  function  method: 


1)  The  initial  design  must  be  in  the  feasible  re¬ 
gion,  i.e.,  C.(p  )  <  0,  j=  1,2,...,  m  and  care  must  be  taken 

J  ~  O  —  “* 

1c 

in  calculating  to  remain  in  the  feasible  region  at  all  times. 
The  augmented  function  loses  its  meaning  outside  the  feasible 
region. 

2)  The  selection  of  an  initial  value  for  the  penalty 
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parameter  r  is  intuitive.  If  the  initial  r  is  too  large,  the 
function  <}>(p,r)  will  be  easy  to  minimize,  but  the  minimum  may 
lie  far  from  the  optimum  point  of  the  original  constrained  pro¬ 
blem.  If,  on  the  other  hand,  the  initial  r  is  too  small,  the 
function  <|>(P,r)  will  be  hard  to  minimize. 

3)  The  minimum  of  the  function  <f>(p,r^)  is  deter¬ 
mined  for  each  r^  value  and  then  r^  is  reduced  to  its  next  value 
r^+^.  This  process  continues  until  the  convergence  criteria 
for  function  F(p)  is  satisfied.  This  nature  of  the  algorithm 
tends  to  increase  the  number  of  iteration  steps  since  smaller 
step  sizes  have  to  be  employed  near  minimum  points  (Ref.  19) . 

As  described  above,  all  three  algorithms  have  cer¬ 
tain  common  features  which  can  be  summarized  as  follows: 

1)  They  all  require  the  calculation  of  the  merit 
function,  the  constraints,  the  merit  function  gradient  and  the 
constraint  gradients  at  every  iteration.  In  the  interior  penalty 
function  approach  this  involves  all  the  constraints;  in  the 
other  two  methods  only  the  currently  active  constraints  are 
involved. 

2)  All  methods  require  a  certain  amount  of  algebra 
in  the  direction-finding  process:  The  calculation  of  the  matrix 
H  for  the  interior  penalty  function  method,  the  linear  pro- 
gramming  solution  for  the  Zoutendijk  algorithm  and  the  calcu¬ 
lation  of  the  matrix  P  for  the  gradient  projection  approach. 

^4 


26 


3)  All  the  algorithms  require  the  solution  of  a  one¬ 
dimensional  minimization  problem,  where  care  must  be  taken  to 
remain  in  the  feasible  region.  For  the  interior  penalty  func¬ 
tion  approach  the  one-dimensional  minimization  involves  the 
calculations  of  the  augmented  function  <f>(p,t)  which  is  always 
non-linear,  whereas  for  the  other  two  methods,  only  the  merit 
function  F(p)  is  involved,  which  can  be  linear.  On  the  other 
hand,  remaining  in  the  feasible  region  requires  very  little 
effort  for  the  interior  penalty  function  method  since  the  con¬ 
straints  are  inherent  in  the  augmented  function,  for  Zoutendijk's 
method  it  involves  an  interpolation  only  when  new  constraints 
are  encountered,  provided  appropriate  0^  values  are  chosen,  for 
the  gradient  projection  method,  assuming  we  are  working  with 
non-linear  convex  constraints,  an  extra  step  is  always  required 
to  come  back  to  the  feasible  region. 

4)  All  methods  involve  certain  parameters  which  can 
only  be  properly  selected  after  enough  experience  with  the 
particular  method  and  with  the  particular  optimization  problem. 
The  penalty  parameters  r^  for  the  interior  penalty  function 
method,  the  scaling  coefficients  0^  for  the  usable-feasible 

ic 

directions  method  and  the  step  size  for  the  gradient  pro¬ 
jection  method. 

The  above  summary  indicates  that  any  choice  of  one  of 
the  algorithms  over  the  other  two  is  going  to  be  problem  depen¬ 
dent.  We  have  the  following  preliminary  data  for  the  aero- 
elastic  optimization  problem: 
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1)  Objective  function  is  linear  in  design  parameters. 

2)  Flutter  and  stress  constraints  are  highly  non¬ 
linear. 

3)  Thickness  constraints  are  linear. 

4)  The  calculation  of  the  flutter  constraint  and  its 
gradients  are  the  major  time  consiming  computational  efforts. 

With  this  much  information  it  is  possible  to  discard 
the  penalty  function  method,  since  it  will  not  take  advantage 
of  the  linearity  of  the  objective  function  and  since  it  will 
include  all  the  constraints  at  all  times.  Also  because  of  its 
nature,  (See  page  25)  it  will  tend  to  involve  a  larger  number 
of  iterations . 

To  make  a  choice  between  the  other  two  is  more  diffi¬ 
cult.  Since  the  objective  function  is  linear,  a  one-dimensional 

* 

minimization  procedure  will  not  yield  the  step  length  a  for 

Si 

either  method.  a  must  be  determined  some  other  way.  This  can 

q 

be  done  empirically  or  can  be  based  on  some  criteria  such  as 

the  reduction  in  merit  function  (Ref.  21) .  In  choosing  a*, 

SI 

the  main  concern  for  Zoutendijk's  approach  will  be  the  violation 
of  new  constraints  whereas  the  main  concern  for  the  gradient 
projection  method  will  be  staying  close  to  the  feasible  region. 
Proper  0j  values  must  be  selected  for  the  Zoutendijk's  approach 
whereas  a  proper  way  of  coming  back  to  the  feasible  region  must 
be  devised  for  the  gradient  projection  method. 

The  final  choice  for  the  optimization  algorithm  in 
this  study  is  the  gradient  projection  concept  which  was  sub- 


28 


sequently  modified  to  handle  the  constraint  violations.  The 
main  reason  for  this  choice  was  to  provide  an  alternative  to 
the  Zoutendijk's  approach  which  has  already  been  successfully 
used  by  several  studies  (Refs.  20,  21)  on  aeroelastic  optimi¬ 
zation  in  recent  years. 


CHAPTER  3 


STRUCTURAL  FORMULATION 

3.1  Sandwich  Idealization 

For  the  examples  considered  here  a  sandwich  type 
structural  idealization  is  used.  This  type  of  structural 
idealization  is  well  suited  especially  for  the  preliminary  de¬ 
sign  of  most  types  of  lifting  surfaces.  This  type  of  ideali¬ 
zation,  as  will  be  seen  later,  is  also  computationally  very 
efficient. 

A  sandwich  plate  (Figure  3.1)  is  the  extension  of  the 
concept  of  I-beam  into  two-dimensions.  It  consists  of  shear¬ 
carrying  core  material  with  negligible  bending  stiffness, 
covered  on  both  sides  with  skins  having  high  inplane  stiffness. 
The  following  assumptions  further  define  the  model: 

1)  t  <<  d  i.e.,  cover  skins  can  be  considered  to  be 
membranes  with  negligible  bending  and  transverse  shear  stiff¬ 
ness. 

2)  A  line  normal  to  the  mid-surface  in  the  unde¬ 
formed  state  remains  straight  and  normal  to  the  mid-surface 
across  the  core  and  the  skins  after  deformation. 

3)  Membrane  stresses  in  the  cover  skins  are  constant 
across  the  thickness  of  the  cover  skins. 

4)  Both  cover  material  and  core  material  are  iso¬ 
tropic. 
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With  these  assumptions,  the  following  equations  ex¬ 
pressing  the  membrane  stresses  in  the  cover  skins  in  terms  of 
second  derivatives  of  the  transverse  displacement  can  be 
written : 
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The  following  equations  give  the  plate  bending 
moments  in  terms  of  second  derivatives  of  the  transverse  dis¬ 
placement: 
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In  both  of  these  equations  v,  the  Poisson's  ratio 
and  E,  the  Young's  modulus,  are  of  cover  skin  material.  From 
equation  (3.2)  we  can  define  the  sandwich  plate  flexural 
rigidity  term  D  as 


D 


E  d2t 
2 (1-v2) 


(3.3) 


In  all  of  the  optimization  problems  discussed  here, 
core  thickness  d  is  held  constant  and  thickness  t  of  both 
upper  and  lower  skins  is  varied.  By  using  a  non-dimensional 
scaling  function  p(x,y)  it  is  possible  to  express  a  new  skin 
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thickness  distribution  in  terms  of  the  initial  uniform  thickness 
distribution  (Figure  3.2); 

t(x,y)  =  p(x,y)tQ  (3.4) 

where  t  is  the  skin  thickness  of  the  initial  design.  From  this 
definition  we  have 

PQ(x,y)  =  1  (3.5) 

Since  the  plate  flexural  rigidity  term  D  is  a  linear 

function  of  t,  (Eq.  3.3)  the  varying  flexural  rigidity  D(x,y) 
can  be  expressed  in  terms  of  the  initial  uniform  flexural 
rigidity  as 

D (x,y)  =  p(x,y)Do  (3.6) 

The  initial  uniform  mass  per  unit  area  of  the  sand¬ 
wich  plate  can  be  divided  into  variable  and  constant  mass  such 
that 

mQ  =  mv  +  mc  =  mQri  +  mQ(l  -  n)  (3.7) 

where  _ 

m 

n  =  —  (3.8) 

m 

o 

Here  variable  mass  is  the  mass  of  the  cover  skins;  constant 
mass  is  the  mass  of  the  core  material.  For  an  actual  wing,  con¬ 
stant  mass  will  include  the  mass  of  the  spars,  ribs  and  fuel. 
Since  variable  mass  is  a  linear  function  of  t,  we  can  write 

m(x,y)  =  mQ[p(x,y)n  +  1  -  hi  (3.9) 

where  m(x,y)  is  the  total  mass  per  unit  area  corresponding  to  a 
new  cover  skin  thickness  distribution. 

As  can  be  seen,  at  any  step  of  the  optimization,  by 
specifying  the  non-dimensional  scaling  function  p(x,y)  both 
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stiffness  and  mass  parameters  of  the  sandwich  plate  are  easily 
defined  in  terms  of  the  initial  stiffness  and  mass  parameters. 

3.2  Finite  Element  Idealization: 

Finite  element  idealizations  are  used  to  discretize 
the  continuous  structural  examples.  The  following  finite  ele¬ 
ments  are  employed: 

1)  Constant  thickness  sandwich  beam  element 

2)  Tapered  sandwich  beam  element 

3)  High-precision  triangular  sandwich  plate  bending 
element  with  linear  thickness  variation. 

Details  of  the  element  characteristics  are  in  Appen¬ 
dices  A,  B,  and  C. 

The  scaling  function  p(x,y)  is  discretized  by  either 
associating  a  with  each  element  (for  constant  thickness  ele¬ 
ments)  or  with  each  nodal  point  (for  elements  with  linear  thick¬ 
ness  variation).  In  this  way  a  continuous  function  p(x,y)  is 
transformed  into  a  design  variable  vector  p.  The  initial  design 
vector  used  in  the  present  studies  will  be 

PoT  =  U,  1,  1,...,  1,  1}  (3.10) 

With  this  approach  it  is  possible  to  write  the  mass 
and  stiffness  matrices  in  the  following  form: 

K(p)  =  l  p.K.  (3.11) 

*  i 

n 

M(p)  =  l  p.M.  +  M 

~  ~  v  x  £  -L  ~  U 

~  X 


(3.12) 
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where  M  is  the  constant  mass  matrix,  K.  and  M.  are  the  indi- 

Z c  ~  1  ~  1 

vidual  contributions  of  the  design  variable  to  stiffness  and 
mass  matrices  respectively.  We  can  further  define  initial  mass 
and  stiffness  matrices  as 

K  =  K  (p  )  (3.13) 

M  =  M (p  )  (3.14) 

7Z  L/  ~  ~  O 

Using  Mq  we  can  write 

M0  =  (1  -  T))M  (3.15) 

which  can  be  used  to  generate  Mc  from  the  initial  mass  matrix. 

The  above  representation  of  the  mass  and  stiffness 
matrices  enables  us  to  generate  mass  and  stiffness  matrices  for 
different  design  variable  vectors  p^.  Element  matrices,  which 
are  independent  of  the  design  variables.  Need  only  be  computed 
and  stored  once  at  the  beginning. 


CHAPTER  4 


AERODYNAMIC  FORMULATION 


This  formulation  basically  follows  Olson's  in  reference 
6.  Two-dimensional,  quasi-steady  aerodynamic  theory  is  used  to 
obtain  the  aerodynamic  forces.  It  is  assumed  that  the  air  flow 
is  parallel  to  the  x-axis  (Figure  4.1)  above  the  panel  and  the 
effect  of  any  air  below  the  panel  may  be  neglected.  In  terms 
of  the  transverse  displacement  function,  W(x,y,t),  the  aero¬ 
dynamic  pressure  acting  on  an  infinitesimal  element  dxdy  is 


Ap  (x,y ,  t) 


2q 
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dxdy 


(4.1) 


where  q  is  the  dynamic  pressure,  Mm  is  the  free-stream  Mach 
number,  and  u  is  the  flow  velocity. 

It  is  possible  to  separate  the  transverse  displace¬ 
ment  function  w(x,y,t)  into  functions  of  space  and  time 


W  (x,y ,  t) 


W (x,y)eyt 


(4.2) 


where,  in  general,  u  is  a  complex  number 

p  =  3  +  iw  (4.3) 

The  borderline  between  stability  and  instability  is  defined  by 
vanishing  3.  At  the  flutter  condition  the  motion  is  given  by 
eiwt  which  is  harmonic. 
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W(x,y,t)  =  W(x,y)eM 


Fig.  4.1  Sandwich  Panel  in  Supersonic  Flow 
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Substitution  of  equation  (4.2)  into  equation  (4.1) 
allows  the  separation  of  the  aerodynamic  force  into 

Ap  =  [Apa  +  Ap(^]e^Jt  (4.4) 


where 


Ap. 


2q _  3w 
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dxdy 


(4.5) 


is  the  primary  aerodynamic  force ,  and 

2 

,  Mr  -2 
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(4.6) 
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is  the  aerodynamic  damping. 

9W  .  .  * 

Ap  is  a  function  of  tt— ;  to  include  its  effect  into 

Si  O  A. 

the  matrix  structural  equation  it  is  necessary  to  derive  an 
aerodynamic  matrix.  This  can  be  done  by  calculating  the  vir¬ 
tual  work: 
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where  w(x,y)  is  the  virtual  displacement.  In  the  finite  ele¬ 
ment  representation  W(x,y)  will  usually  have  the  following  form: 


W(x,y>  =  fT  Te  «e 


(4.8) 


where  f  is  the  vector  of  interpolation  functions  in  x  end  y  j 

T  is  the  transformation  matrix  that  relates  interpolation  func- 
-e 

tion  coefficients  to  the  element  generalized  displacements ,  and 

VJ  j_s  the  vector  of  element  generalized  displacements r  which  are 
~e 

element  nodal  point  displacements. 
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Introducing  equation  (4.8)  into  equation  (4.7)  gives 
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Here  Ae  is  the  consistent  aerodynamic  matrix  for  an  element. 

Ap^  is  a  function  of  W(x,y).  Its  effect  can  be  in¬ 
cluded  by  deriving  an  Aerodynamic  Damping  matrix,  De*  Follow¬ 
ing  the  same  steps  as  above,  we  can  obtain 
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It  should  be  noted  that  the  double  integral  in  equation 
(4.13)  is  the  same  as  the  one  required  for  the  derivation  of  the 
element  mass  matrix  for  constant  mass. 


CHAPTER  5 


CONSTRAINTS 

In  this  chapter  equations  expressing  the  flutter, 
stress  and  thickness  constraints  in  terms  of  the  design  variable 
vector  p  are  derived. 


5.1  Flutter  Constraint 


In  order  to  express  the  flutter  constraint  in  terms  of 
p,  it  is  necessary  to  write  the  equation  of  motion  for  the  sys¬ 
tem.  In  the  finite  element  approach  this  can  be  derived  from 
the  equation  of  motion  for  an  element. 

The  element  stiffness  can  be  represented  by  the  ele¬ 
ment  stiffness  matrix  Ke  which  can  be  derived  from  strain  energy. 
It  will  be  a  function  of  p(x,y)  which,  in  turn,  can  be  expressed 
in  terms  of  the  elements  of  p. 

The  contribution,  AI,  of  the  inertia  force  (Figure 


4.1)  , 


2  - 

AI  =  mQ(pri  +  1  -  n)u  W  dxdy 


(5.1) 


can  be  represented  by  the  element  mass  matrix,  M  ,  which  can  be 
derived  from  either  kinetic  energy  or  from  virtual  work.  It 
will  also  be  a  function  of  p(x,y).  References  15  and  16  can  be 
mentioned  as  good  sources  for  detailed  derivations  of  mass  and 
stiffness  matrices. 

The  contributions  of  Ap  ,  the  primary  aerodynamic 
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force  (Equation  4.5),  and  Ap^,  the  aerodynamic  damping  (Equa¬ 
tion  4.6),  can  be  represented  by  the  element  aerodynamic  matrix, 
A£,  and  element  damping  matrix,  Dg,  derivations  of  which  have 
already  been  discussed  in  Chapter  4. 

We  can  write  the  element  equation  of  motion  in  terms 
of  these  matrices; 


Fe  =  [Ke(p)  +  U  Me(p)  +  Ae  +  u  De]we 


(5.2) 


where  Fe  and  W£  are  generalized  forces  and  displacements 
respectively.  They  are  related  to  the  element  nodal  point  de¬ 
grees  of  freedom.  We  can  rewrite  equation  (5.2)  by  using  non- 

dimensional  matrices  K  ,  M  ,  A  and  D  ; 

~e  ~e  ~e  ~e 
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where  £  is  a  convenient  dimension  used  in  scaling  all  finite 

element  dimensions  when  deriving  K  .  M  ,  A  and  D  .  It  should 

-e  ~  “ 

be  noted  that  5  =  M  for  p  =  1.  The  matrix  coefficients  in 

«e  «e 

equation  (5.3)  can  be  rearranged  such  that  we  have 

?e  =  Do[?e!p)  +  <S“)2  Se(p)  +  “ole  +  9ao  O  °e)>?e  (5.4) 
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is  a  convenient  frequency  scale  factor 


(m;  -d1/2  °o 


is  the  non-dimensional  dynamic  pressure  parameter,  and 


=  2a 


M2  -  2  , 

°°  1 

(M2  -  1 )3//2  mo“o 

00 


(5.6) 


(5.7) 


is  the  non-dimensional  aerodynamic  damping  parameter. 

The  matrices  K  ,  M  ,  A  and  5  .  and  vectors  F  and 
~  0  ~  e  %  6 

W  can  be  assembled  in  the  usual  manner  to  form  the  system  equa- 

~e 

tions  of  motion.  For  no  external  forces,  these  equations  can 
be  written  as 


[K  (p)  +  (Jp)2  M(p)  +  aoA  +  gao  (Jf*)  D]W  =  0 


(5.8) 


where  W  is  the  vector  of  system  generalized  displacements,  which 
are  related  to  the  system  nodal  point  degrees  of  freedom.  Here 
K(p)  and  M(p)  have  the  forms  expressed  in  equations  (3.11)  and 
(3.12).  It  should  be  also  noted  that 


D  =  M 
~o 


(5.9) 


where  M  is  the  same  as  in  equation  (3.14). 

~o 

For  q  j-  0  Equation  (5.8)  can  be  recognized  as  a 
^ao 

quadratic  eigenvalue  problem  in  (~p) •  We  can  reduce  this  equa- 

o 

tion  into  a  set  of  coupled  first  order  equations  by  introducing 


X  =  y/w. 


(5.10) 


V  =  xw 


(5.11) 
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and 


U  = 


(5.12) 


Using  equations  (5.10)  and  (5.11)  in  equation  (5.8) 

we  get 

[K  +  a  A]  W+gDV  =  -  X  M  V  (5.13) 

and  from  equation  (5.11)  we  have 

Y  =  A  W  (5.14) 


or 


-M  [K  +  ot  A]  W  - 

~  Z  O  ~ 

g  M-1 

D  V 

=  x  y 

(5.15) 

and 

V  = 

X  W 

(5.16) 

Equations  (5.15) 

and  ( 5 . 

16) 

can  be  rewritten  by  using 

U  of 

equation  (5.12); 

C  U 

=  X  U 

(5.17) 

where 

^ao 

M-1  D 

1 

1 

1 

1 

-  M_1[K  +  a  A] 

*vr  jv  O  r* 

(5.18) 

c  = 

I 

I - 

1 

0 

1 

I 

Equation  (5.17)  is  a  linear  eigenvalue  problem  of  size 
2N  x  2N  where  N  is  the  number  of  system  degrees  of  freedom.  For 
under  damped  systems  (see  ref.  24  for  critical  damping,  under 
and  over  damped  systems)  equation  (5.17)  has  N  complex  conjugate 
eigenvalues  and  corresponding  N  complex  conjugate  eigenvectors. 


43 


Recall  that  from  (5.12),  the  last  half  of  the  eigenvectors 
(i.e.,  W)  are  modes  shapes  of  motion,  while  the  eigenvalues  are 
the  corresponding  complex  frequencies  scaled  by  uq  which  is 
real  and  positive.  In  the  complex  frequency  the  real  part  is 
the  damping  part,  the  imaginary  part  is  the  frequency  of 
oscillatory  motion  (equations  4.2  and  4.3).  Complex  mode  shapes 
indicate  that  there  are  phase  differences  between  different 
points  of  the  structure  in  motion. 

In  Chapter  4  it  has  been  pointed  out  that  the  border¬ 
line  between  stability  and  instability  is  defined  by  the  real 
part  of  y  becoming  zero.  Since  X  =  y/a)Q,  the  eigenvalues  with 
positive  real  parts  indicate  aerodynamic  instability,  or  flutter. 
This  gives  us  the  flutter  constraint  as; 

Cf  =  Real  (X)  <  0  (5.19) 

where  X  is  any  eigenvalue  of  the  equation  (5.17). 

For  a  given  damping  parameter  gaQ  (it  should  be  pointed 
out  that  any  existing  structural  damping  can  be  represented  by 
adding  a  structural  damping  parameter  to  the  aerodynamic  damp¬ 
ing  parameter) ,  and  matrices  K,  M,  D  and  A,  the  flutter  con- 

z  ~  z  ~ 

dition  can  be  reached  by  increasing  aQ  from  zero.  The  value  of 

a  that  causes  flutter  is  called  a  .  Since  the  final  minimum- 
o  cr 

weight  design  is  required  to  have  the  same  flutter  speed  as  the 
initial  design  we  have  the  matrices  K  ,  M  #  D  and  A  and  the 
relationship 

MED 
~o  ~ 
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which  reduces  the  C  matrix  to; 


So 


o 


(5.20) 


For  aQ  =  0,  the  eigenvalues  of  the  CQ  matrix  are  sim¬ 
ply  the  damped  natural  frequencies  of  the  system.  The  real 
parts  of  all  the  eigenvalues  will  be  “9ao/2*  As  ao  ^ s  increased; 
magnitudes  of  the  imaginary  parts  of  two  complex  conjugate  pairs 
get  closer  (generally  the  lowest  two  frequency  pairs  are  in¬ 
volved)  .  The  real  parts  remain  the  same  until  the  imaginary 
parts  become  nearly  equal  (a  required  condition  for  energy  trans¬ 
fer  from  one  mode  to  the  other).  After  that,  the  imaginary 
parts  reamin  equal  but  increase  in  magnitude,  the  damping  of  one 
pair  increases  as  the  damping  of  other  pair  decreases.  At  aQ 
=  a  one  pair  will  have  g  ,  the  other  pair  will  have  zero 
damping,  which  means  the  motion  at  flutter  is  pure  harmonic. 

For  aQ  >  acr  the  same  trend  continues;  the  real  part  of  the 
flutter  pair  becomes  positive  (indicating  instability) ,  the 
other  pair  gets  further  damped,  such  that  the  equation 


Real  (A^  +  Real  U2)  =  ~gao  (5.21) 

is  always  valid  (see  figure  5.1). 

As  a  is  increased  from  zero  to  a _  the  real  parts 

o  cr 

of  all  the  other  pairs  remain  at  ~9ao/2'  the  imaginarY  parts 
change  slightly. 
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The  numerical  determination  of  acr  for  the  initial 
uniform  design  requires  solving  equation  (5.17)  for  different 
a  values  with  C  matrix  of  equation  (5.20).  Only  the  eigen- 
values  need  to  be  calculated.  By  looking  at  the  real  parts  of 
the  eigenvalues  it  is  possible  to  close  in  on  acr  using  a  bisec¬ 
tion  method. 

During  the  optimization,  in  order  to  calculate  the 
flutter  constraint,  matrix  C  of  equation  (5.18)  has  to  be  formed. 
The  aerodynamic  matrix  A  is  already  generated  for  the  preliminary 
flutter  analysis,  the  damping  matrix  D  is  obtained  by  storing 
Mq.  Matrices  M  and  K  need  to  be  generated  at  every  design  step 
according  to  the  equations  (3.11),  (3.12)  and  (3.15).  In  the 
problems  considered  here,  the  flutter  constraint  is  kept  active 
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at  all  design  steps.  The  eigenvalues  of  the  C  matrix  are 
calculated  by  first  reducing  the  matrix  into  a  Hessenberg  form 
(Ref.  18)  and  then  applying  QR  iterations  (Ref.  18)  on  the 
Hessenberg  form.  Among  the  eigenvalues  the  pair  with  the  lar¬ 
gest  positive  real  part  is  chosen  to  be  the  flutter  eigenvalue. 
For  cases  where  all  the  eigenvalues  have  negative  real  parts , 
the  pair  with  the  smallest  imaginary  part  in  magnitude  is  chosen 
to  be  the  flutter  eigenvalue. 

Other  points  of  interest:  a)  During  the  optimization 
process,  this  way  of  expressing  the  flutter  constraint  enables 
us  to  detect  flutter  involving  modes  other  than  the  original  two 
modes,  b)  The  flutter  constraint  is  highly  non-linear  in  p;  the 

functional  relationship  of  the  matrix  C  to  p  is  through  matrices 

*  ~ 

K  and  M  which  are  functions  of  p  as  expressed  in  equations 
(3.11)  and  (3.12). 


5.2  Stress  Constraint 

Stress  constraints  are  considered  only  for  the  speci¬ 
fied  uniform  static  loading  conditions  used  as  a  first  approxi¬ 
mation  to  the  aeroelastic  loading. 

Von  Mises  yield  criterion  is  used  to  obtain  stress 
constraint  equations.  For  cover  skin  stresses  where 

a  =  a  =  a  =  0,  the  Von  Mises  yield  criterion  reduces  to 
zz  zy  zx  '  J 


i  [  (a  -  a)2  +  a2  +  CT2  +  6r2]=K2 
6  x  y  y  x  xy 


(5.22) 
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Here  K  is  yield  strength  in  pure  shear  or 

1 


K  = 


V  3  ^ayield  pure  tension)  (5.24) 


Using  equation  (3.1)  for  stresses  we  can  write  equa¬ 
tion  (5.23)  in  terms  of  curvatures  W  ,  W  and  W  ; 

xx  yy  xy 


Pr?  2 

- 2[ {Wxx  +  vWyy)  "  (Wxx  +  vWyy}  (vWxx  +  V 


2(1  -  v) 


+  (vWxx  +  Wyy)2  +  3  (1-v) 2  Wxy2]  =  3K2 


(5.25) 


We  can  reduce  equation  (5.25)  to  a  standard  constraint 


form: 


(v2  -  V+l) (W2X  +  Wyy)  -  (v2 


4v  +l)WxxWyy 


+  3(1  -  v) 2W2y 


A  <  0 


(5.26) 


where 


6K2 (1-v)2 
Ed 


(5.27) 


In  the  finite  element  analysis,  for  any  point  on  the 

structure  where  a  stress  constraint  is  applied,  curvatures  Wxx, 

W  and  W  can  be  calculated  from  the  element  nodal  point  dis- 
yy  xy 

placement  vector  W£.  Using  equation  (4.8)  we  can  write 


1 

w 

[fT  1 

XX 

-XX 

W 

•  _ 

fT 

yy 

~yy 

w 

f  l. 

xy  J 

-xy 

(x,y) 

—  - 

TW 

»e~e 


(x,y) 


(5.28) 
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Wg  can  be  obtained  from  the  system  nodal  point  dis¬ 
placement  vector  W,  which  in  turn  can  be  calculated  from  the 
equation  of  equilibrium  for  static  loading; 

K(p)W  =  L  (5.29) 


where  K  is  the  system  stiffness  matrix  and  L  is  the  load  vector, 
From  equations  (5.28)  and  (5.29)  we  get 


f  \ 

W„ 

fTl 

XX 

-XX 

W 

fT 

yy 

' 

~yy 

w 

fT 

xy 

-xy 

k  J 

(x,y) 

L  — 

T  X  W 


(5.30) 


(x,y) 


where  the  rectangular  matrix  X  picks  up  the  vector  W  from  W 


i.e: 


v?  =  x^  w 

~e  ~e  ~ 


(5.31) 


It  should  be  noted  that,  in  equation  (5.30)  only  W 
is  a  function  of  p  (according  to  equation  (5.29)).  This  fact 
will  be  used  in  deriving  the  equations  for  stress  gradients  in 
Chapter  7 . 

Stress  constraints  are  only  considered  for  the  delta 
wing  problem.  The  load  vector  is  calculated  for  a  uniform  load 
distribution  by  condensing  loads  at  the  nodal  points.  The  in¬ 
tensity  of  the  uniform  load  distribution  is  determined  to  pro¬ 
duce  yield  for  the  uniform  initial  design.  The  constraints  are 
applied  at  the  nodal  points  where  transverse  displacement 
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second  derivatives  are  obtained  directly  from  nodal  point  de¬ 
grees  of  freedom  (see  Chapter  9  for  details) . 

5.3  Thickness  Constraint 

Thickness  constraints  are  applied  on  the  elements  of 
the  design  variable  vector  p.  Since  negative  thickness  is  not 
physically  possible,  the  condition  pi  >_  0  must  be  satisfied  by 
all  p^  at  all  stages  of  the  optimization.  The  following  is  the 
standard  constraint  equation  for  the  ith  design  variable 

C.  .  =  p_.  “  P„.  <  0?  0  <  p  .  <1  (5.32) 

For  cases  where  nodal  point  cover  skin  thickness  para¬ 
meters  used  as  design  variables,  the  above  equation  is  sufficient 
to  constrain  thickness  within  elements  since  linear  thickness 
distribution  is  assumed. 

pmin  =  0*1  is  used  in  the  optimization  examples. 


CHAPTER  6 


GRADIENT  PROJECTION  ALGORITHM  WITH 
MULTIPLE  CONSTRAINTS 


The  basic  ideas  of  the  gradient  projection  method  have 
been  discussed  in  Section  2.4  in  forming  a  comparison  with  the 
other  popular  algorithms.  It  was  pointed  out  that  the  main  pro¬ 
blem  with  the  algorithm  was  to  stay  on  or  within  the  convex  con¬ 
straints.  The  equations  for  the  algorithm  used  in  this  study, 
where  the  minimization  of  the  objective  function  and  returning 
to  violated  constraints  are  performed  in  the  same  step,  will  now 
be  derived. 

We  assume  that  at  the  design  point  p  there  are  m 

c 

critical  constraints,  which  for  convenience  are  considered  to 
be  numbered  1,  2,  ...,  m  ,  i.e.  constraints  for  which 


C.  (p  )  >  0  ,  j  =  l,  ...,in 

j  '  —  9  c 


(6.1) 


we  want  to  calculate  a  new  design  point  Pg+^  where; 


F(eq+i>  <  F<eq) 


(6.2) 


and 


Cj (6q+l)  -  0  '  3  =  ,  mc 


(6.3) 


We  define  Ap^  as  the  vector  required  to  take  us  from 


fig  t0  ~q+l '  1-e‘ 
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eq+i  =  eq  +  Aeq  <6-41 

Using  linear  approximations  for  F(p)  and  Cj (p) , 

j  =  1,  m  at  p  we  can  write  AF,  the  change  in  the  objec- 

c  f 

tive  function,  and  AC.,  j  =  1,  . ..,  m  ,  changes  in  the  critical 

J  C 

constraints  along  Ap  as; 

~  Si 

AF  =  Ap  T  VF (p  )  (6.5) 

~q  ~q 

AC.  =  Ap  T  VC. (p  ) ,  j  =  1,  ...,  m  (6.6) 

J  ~q  j  ~q 

We  define  a  norm  for  Ap  ; 

Rq2  =  46qT  A2q  (6-7) 

where  R  ,  the  step  size,  is  a  pre-selected  positive  number  for 
this  particular  design  step  (see  Sections  8.3,  9.4,  9.5  and 
Chapter  10  for  the  selection  of  the  step  size) . 

Now  the  calculation  of  Ap  can  be  cast  as  an  optimi- 
zation  problem  with  equality  constraints: 

Maximize  (-AF) 

Subject  to  AC.  =  -  C. (p  ) ,  j  =  1,  ...,  m 

j  j  ~q 

and  Rg2  -  AfiqT  Afiq  =  0 

so  that  the  requirements  in  equations  (6.2)  and  (6.3)  are  met  to 
within  a  linear  approximation. 

Using  Lagrange  multipliers  y  j ,  j  =  1/  •••»  mc  an<^  vr 
we  define  the  following  functional: 


52 


m 


r  -  -AeqT  VF<eq>  +  |  Vi  [AecrT  vcj  <e<j)  +  <Ve„>i 


d  ~q 


+  V-fR  2  -  Ap  T  Ap  ) 
R  q  t-q  S-q' 


sr 


0  gives 


m 


-VF(eq>  +  l  VW  '  2Veq  -  ? 

ApQT  VCj(pq)  +  C^(po)  =  0,  j  =  1,  ..., 


j  ~q 


m 


(6.8) 


(6.9) 

(6.10) 


2  T 

R  Ap  Ap  =0 

q  ~q  ~q 


(6.11) 


We  can  rewrite  equations  (6.9)  and  (6.10)  using  matrix 
notations  as  follows: 


-VF(eq>  +  gq  H  -  2',RAPg  =  ? 


(6.12) 


where 


G 

~q 


c 

~q 


G  Ap  +  C  =0 
»q  ~q  ~  q  ~ 


*  [VCl<2q>'  VC2(Bq) . VCmc(Sq>] 

tylf  y2 '  **•'  ymc^ 

*Cl*2q*'  C2  ^6q^  '  Cmc(2q)] 


(6.13) 

(6.14) 

(6.15) 

(6.16) 


From  equation  (6.12)  we  obtain 


% 


2v 


R 


[-VF  (p  )  +  G  y] 


-q 


zq  ~ 


(6.17) 
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which  can  be  substituted  in  equation  (6.13)  to  give 


2^gqT'-w'eg>  +  fqHi  +  ?q  -  2 


(6.18) 


G„T  G  y  =  -2vp  C  +  G  T  VF (p  ) 


(6.19) 


Assuming  columns  of  G  are  linearly  independent,  the 
T 

matrix  [G  G  ]  is  nonsingular  and  can  be  inverted.  Thus  equa- 
tion  (6.19)  yields: 

y  =  -2vp [G  ^  G]  C  +  [G  T  G]  GT  VF (p  )  (6.20) 

~  R  »q  ~q  ~q  »q  «q  «q  ~q 


Using  this  expression  for  y  in  equation  (6.17)  we  get 

Ap  =  (G  [G  T  G  ]  G  T  -  I]  VF  (p  ) 

~q  2Vp  l~ql"q  ~qJ  ^q  »J  vKqy 


-G  [G  G  ]  C 
~q  «q  »qJ  ~q 


(6.21) 


vR  can  be  calculated  by  using  equation  (6.11).  For  simplicity 
we  can  write 


e  +  f 


(6.22) 


where  e  =  i  [Grt[G„T  G  ]  GaT  -  I]VF(p  ) 
2  ~q  x:q  ~q  ~q  »  ~q 


(6.23) 


!  =  -gqtGq  G  1  Cg 


(6.24) 


Now  equation  (6.11)  can  be  written  as 


(i-  eT  +  fT)  (~  e  +  f) 
VR  ~  VR  ~ 


0 


(6.25) 
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which  reduces  to 


where 


2  2 

‘V  '  clvR 


a  =  e  e 


b  =  2e  f 


c  =  f  f 


-  bvR  -  a  =  0 


(6.26) 

(6.27) 

(6.28) 
(6.29) 


Using  expressions  for  e  and  f  we  get 


b  =  VF(p  )T[G  [G  TG  ]  1  G  T  -  I]  [-Ga[GaTGo] 

4  ^4  s;4  ~4  ^4  s  ^4  ~q 


(6.30) 


or 


b  =  VF (p  )  [0]C^  =  0 

~4  ~  ~4 


(6.31) 


which  means  that  e  is  orthogonal  to  f.  Using  this  information 
in  equation  (6.26)  we  get 


VR  (  ^  2 


) 


1/2 


Rq  "c 


and  finally  Ap^  becomes 


R  -  c  1/2 

Ap„  "  ^ — a - >  e  +  f 

a.  ~  ~ 


(6.32) 


(6.33) 


where  all  the  terms  have  already  been  defined. 
It  should  be  noted  that 


2e  =  S 

~q 


(6.34) 


where  S  has  been  defined  in  equation  (2.13).  So  we  have  the 
design  change  vector  Ap  as  being  composed  of  two  orthogonal 
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R  -c  1/2 

vectors,  e'  =  (-3 - )  e,  and  f,  where  e’  is  in  the  projected 

gradient  direction  along  which  the  changes  in  constraints  are 
zero  and  f  is  in  an  orthogonal  direction  along  which  the  changes 
in  constraints  are  equal  to  the  constraint  violations  at  in 
magnitude  with  opposite  signs,  that  is 


(6.35) 


and 


(6.36) 


For  cases  where  the  norm  of  the  vector  f  (i.e.  c)  is 
bigger  than  R  2,  the  equation  (6.33)  cannot  be  used  for  the  cal- 
culation  of  Ap  ,  since  the  term  inside  the  square  root  becomes 

~q 

negative.  For  such  cases  the  following  relationship  is  used  to 

calculate  Ap  ; 

~  Vi 


A  f 

/c~  ~ 


(6.37)' 


where  only  the  component  required  to  come  back  to  the  violated 
constraints  is  employed  and  the  relationship 

Eq  =  Aeq  % 


i 


is  still  valid. 

Figures  6.1  and  6.2  illustrate  the  design  steps  given 
by  equations  6.33  and  6.37  on  a  two-dimensional  design  space. 

For  the  example  problems  presented  in  Chapters  8  and 


Merit  Function 


Feasible 

Region 


f  AF(Pq> 


C1  =  0 


7<W 


*2q  = 


R  -  c  1/2  \  e 

(-9 - )  e  + 


c  <  Rq 

Arc  Radius  =  R„ 


Fig.  6.1  Primary  Design  Change  Vector 


Merit  Function 


<^ci  -  0 


Feasible  Region 


Kq+l 


'  VF  ( P  ) 
~q 


c  >  R 


Arc  Radius  =  F 


VC1  ^q* 

Fig.  6.2  Secondary  Design  Change  Vector 
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9,  the  flutter  constraint  is  kept  active  at  all  design  steps. 
Since  aQ  is  determined  to  produce  flutter  for  the  initial  uni¬ 
form  design,  pQ  is  on  the  flutter  constraint.  For  design  steps 
where  <  0,  is  set  equal  to  zero  when  constructing  the 

vector  C  . 

~q 

For  the  delta-wing  problem  the  stress  constraints  are 

considered  to  become  active  when  the  condition  C  >  -  e  is 

s 

satisfied  for  any  nodal  point.  Here  e  is  a  positive  parameter 
which  is  made  successively  smaller  as  the  optimum  design  is 
approached  (see  section  9.4).  For  cases  where  -e  £  Cg  £  0,  Cg 
is  set  equal  to  zero  when  constructing  the  vector  C  . 

Thickness  constraints  become  active  when  there  is  a 
violation  (i.e.  pi  £  Pmin) • 

It  should  be  noted  that  the  derivation  of  Ap^  as  out¬ 
lined  above  does  not  include  any  provision  for  stress  or  thick¬ 
ness  constraints  which  are  not  active  for  p  but  becomes  active 
for  Pq+1*  This  fact  created  some  problems  in  the  earlier  runs 
for  the  delta  wing  optimization;  mainly  a  stress  or  thickness 
constraint  which  was  not  active  at  a  particular  p^  became 
severely  violated  at  Pq+1  and  it  took  several  iterations  to 
bring  the  design  back  into  the  feasible  region  during  which  the 
flutter  constraint  had  to  be  calculated  several  times.  Since  the 
calculation  of  the  flutter  constraint  is  the  main  time-consum¬ 
ing  computational  effort,  the  following  procedure  is  used  to 
reduce  the  computation  time:  After  calculating  Pg+^  the  stress 
and  thickness  constraints  were  checked,  if  there  was  a  new  stress 
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or  thickness  constraint  violation,  a  parameter  aD(0  <  a_  <  1)  is 

calculated  using  linear  interpolation  to  scale  Ap  so  that 

~q 


^q+1 


=  eq  +  “r  Aeq 


(6.38) 


is  in  the  feasible  region  (see  Figure  6.3  for  a  2~dimensional 
illustration).  The  procedure  proved  to  be  very  successful. 

Other  points  of  interest  about  the  algorithm  and  its 
application  to  the  aeroelastic-stress  optimization  are: 

1)  The  method  maximizes  (-AF)  subject  to  returning 

back  to  the  feasible  region,  but  there  is  no  guarantee  that  the 

condition  F(Pq+^^  <  F(Pg)  wiH  be  satisfied,  especially  for 

cases  where  Ap  is  required  primarily  for  coming  back  to  con- 
~  SI 

straints.  However  the  application  of  the  algorithm  showed  that 
this  fact  does  not  create  problems  in  practice. 

2)  Considering  the  fact  that  the  structural  weight  is 
a  linear  function  of  the  design  variables,  equation  (6.5)  is  an 
exact  expression,  however  the  equations  (6.6)  are  exact  only  for 
linear  constraints  (e.g.  thickness  constraints)  and  are  approxi¬ 
mate  for  non-linear  constraints  (e.g.  flutter  and  stress  con¬ 
straints)  .  Their  accuracy  for  non-linear  constraints  is  a 

function  of  R  ,  the  step  size,  and  the  curvature  of  the  con- 
SI 

straints . 

3)  Although  the  algorithm  is  capable  of  handling  con¬ 
straint  violations,  getting  far  away  from  the  feasible  region 
would  be  undesirable  since  the  gradient  information  obtained  at 
such  points  would  be  misleading. 


6.3  Reduction  in  the  Design  Change  Vector  Due  to  a 
New  Constraint  Violation 
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4)  From  2)  and  3)  it  becomes  apparent  that  the  choice 

of  R  is  extremely  important.  A  very  large  value  for  R  could 
4  q 

make  the  algorithm  unstable.  On  the  other  hand  a  very  small 

value  for  R  would  mean  an  increased  number  of  iterations  for 
Si 

convergence. 

5)  The  change  in  the  merit  function  from  one  design 
to  the  next,  is  used  as  the  convergence  criteria  for  the  appli¬ 
cations  (see  results  in  Chapters  8  and  9) . 

6)  For  a  single  constraint  the  above  equations  can  be 
simplified  by  replacing  with  VC(pg).  Equation  (6.21)  be¬ 
comes 


2^  i|  7c(pq)  -  w<eq>)  -  Is  vc(eq> 


(6.39) 


where 


A  =  VC(p  r  VF(p  ) 
^  Si  ^  SI 


B  =  VC (Pg) ^  VC(pq) 


can  be  calculated  using  equation  (6.32).  Employing 


equation  (6.39)  we  now  have 


1  r_  A  i 

a  ~  J  1°  ~  B~] 


(6.40) 


c  = 


(6.41) 


where 


D  =  VF(p  r  VF(p  ) 

~  ~  Si 
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Fig.  6.4  Flow  Chart  for  the  Gradient  Projection  Algorithm 


62 


These  equations  yield 


and  equation  (6.39)  becomes 


BR  2  -  cl  1/2  c 

Ap  =  (—a - gi)  [g  VC(p  )-  VF(p  )]  -  vc(p)  (6.43) 

^  BD  -  A  Vi  ~4  °  ~q 


Figure  6.4  shows  a  simplified  flow  chart  of  the 
algorithm  based  on  the  above  developments . 


CHAPTER  7 


GRADIENT  CALCULATIONS 

The  gradient  projection  algorithm  developed  in  Chap¬ 
ter  6  requires  the  calculation  of  the  gradients  of  the  objective 
function  and  the  gradients  of  the  critical  constraints.  In 
this  chapter  the  analytical  expressions  used  for  these  calcu¬ 
lations  are  derived.  ! 

7.1  Gradients  of  the  Objective  Function 

It  has  already  been  indicated  that  the  objective 
function  used  in  this  study  is  linear  in  the  design  variables. 
This  means  that  the  gradient  vector  VF(p  )  is  a  constant  vector 
everywhere  in  the  design  space.  It  needs  to  be  calculated  only 
once  and  stored.  The  expressions  for  objective  functions  and 
the  gradients  are  given  in  the  text  whenever  the  need  arises 
(Chapters  8  and  9) . 

7.2  Gradients  of  the  Flutter  Constraint 

The  flutter  constraint  was  given  in  equation  (5.19)  as 

Cf  =  Real  (X)  £  0 

where  X  was  determined  by  the  complex  eigenvalue  problem  (5.17) . 

In  equation  (5.17)  the  coefficient  matrix  C  contains  the  mass 

and  stiffness  matrices  M  and  K,  which  are  functions  of  p.  For 

3 

the  gradient  calculations  we  need  the  terms  [Real  (X)] 
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expressed  analytically. 

Equations  derived  by  Rogers  (Ref.  11,  similar  equa¬ 
tions  were  given  by  Van  de  Vooren  in  Ref.  12)  are  employed  to 
achieve  this.  For  the  general  eigenvalue  problems  with  real  or 
complex  eigenvalues  and  eigenvectors 

ct^A(x)(jK  +  B(x)<(k  =  0  (7.1) 


and 


a.AT(x)6.  +  BT(x)9.  =  0 

1  ~  ~  1  ~  ~  1  ~ 


(7.2) 


where  the  right  and  left  hand  eigenvectors  are  normalized  such 
that: 


0  .  A  <f> .  =  1 


(7.3) 


Rogers  (11)  gives  the  following  expression 


ai,X  *  -  ?i  l“iJ,X  +  S,x'h 


(7.4) 


where  ,x  means  differentiation  with  respect  to  x. 

The  flutter  equation  (5.17)  can  be  put  into  the  forms 
of  equations  (7.1)  and  (7.2).  We  start  from  equations  (5.13) 
and  (5.14)  and  rewrite  them  as: 


g  DV  +  [K  +  a  A]W  =  -  X  M  V 

yao«~  «  O  ~  ~  a  ~ 


-  M  V  =  -  XMW 


which  can  be  reduced  to 


M  j  0 

1  <v 

| 

9ao  B  1 

1 

K  +  aA 
^  o  ~ 

1 

u  + 

- - 4  — - — - 

6"|  M 

-  I  - 

-M  , 

-  i 

L  i 

0 

U  =  0 


(7.5) 

(7.6) 


(7.7) 
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where 


U  = 


y 

W 


as  defined  in  equation  (5.12). 

Equation  (7.7)  is  of  the  same  form  as  equation  (7.1). 
By  transposing  equation  (7.7)  we  get 


M 

\ 

1 

1 

r~' 

0 

u  + 

—  1 

g  D  J 

ao  »  ! 

rn  "1 

-M 

Q 

i 

i 

<  ?? 

O 

3 

+ 

_ 1 

0 

which  is  analogous  to  equation  (7.2).  Here  the  symmetry  pro¬ 
perty  of  matrices  M,  D  and  K  are  used.  A  is  not  symmetrical. 

*  z  »  * 

Equation  (7.8)  can  be  reduced  to  the  form 


C  U 


X  U 


(7.9) 


where 


-M_1 [K+a  AT]  !  0 

*  «  o»  i  ~ 


(7.10) 


from  which  the  left  hand  eigenvectors,  U,  can  be  calculated 
using  the  eigenvalues  of  C  matrix. 

Similar  to  equation  (7.3)  we  can  write 


M 

1 

| 

0 

-T 

1 

1 

5; 

U 

— 

“1 - 

_  - - 

0 

1 

1 

M 

=  1 


(7.11) 


by  which  the  left  and  right  hand  eigenvectors 
Using  the  analogy  of  equation  (7.4) 
that  only  M  and  K  are  functions  of  p,  we  get 


can  be  normalized, 
and  remembering 
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X. 

ifp- 


-T 

U. 


M 


'Pj 


0 


0 


M 


,pj 


o  j 

K 

Z  I 

*fp 

■t  - 

— 

-M  ! 

0 

L-pj  ! 

U. 


(7.12) 


substituting 


Ai  =  +  1  a»i)/(0o  =  ^  +  i  wi  (7.13) 


reduces  equation  (7.12)  to 


3.  +  ico.  =  -  UT 

1/ Pj  ifPj  ~i 


3  -M  I  K 

, p •  j  ~ , p  . 
- J_i - J 

-M  !  3,-M 
~  f  P  j  |  1  ^  f  P  j 


+  i 


a) .  M  !  O 

l~fpj  * 

O  !  w  .  M 

*  l  l«/P 


Hi 


1 

(7.14) 


When  right  hand  and  left  hand  eigenvectors  U  and  U  are 

calculated  from  equations  (5.17)  and  (7.9)  and  normalized 

according  to  equation  (7.11),  the  real  part  of  the  equation 

(7.14)  will  give  us  the  required  gradients. 

Using  equations  (3.11)  and  (3.12)  for  K  and  M,  we  can 

*>  ~ 

write: 


K 


(7.15) 


and 


M  =  M.  (7.16) 

« /  P  j  ~3 


where  K.  and  M.  are  the  individual  contributions  of  the  design 
~  J  ~  J 
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variable  to  the  stiffness  and  mass  matrices  respectively. 

These  matrices  can  be  assembled  from  element  mass  and  stiffness 

matrices  in  a  process  similar  to  assembling  system  mass  and 

stiffness  matrices.  This  time,  however,  only  the  elements 

associated  with  a  particular  will  contribute. 

-  T 

It  should  be  noted  that  C  is  not  C  ,  since 


-g  M  1  D 

aO  ~  s: 


-[K  +  a  AT]M-1 

%  o  ~  ~ 


(7.17) 


which  is  different.  Even  though  C  will  have  the  same  eigen¬ 
values  as  C,  its  eigenvectors  are  different  from  those  of  C. 

In  practice  the  right  and  left  hand  eigenvectors  are 
calculated  by  employing  inverse  power  iteration  (Ref.  18)  using 
flutter  eigenvalue  as  the  spectrum  shift. 


7.3  Gradients  of  the  Stress  Constraints 


yields 


Differentiation  of  equation  (5.26)  with  respect  to  p^ 


Cs,e.  -  2<"  -  »  +  ♦VWJ> 

-  (V2  -  4v  +  lXW^W  +  WyyWXX(p.) 


+  6 (1  -  V)  W  W 

xy  xy,p. 


(7.18) 


The  gradients  of  the  curvatures  can  be  obtained  from 


the  gradients  of  the  element  nodal  point  displacement  vector; 
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similar  to  equation  5.28  we  can  write 


xx,  p. 


yy  r  p. 


fT 

~Yy  *e~e  ,p. 


(7.19) 


xy,  p. 


(x,y) 


(x,y) 


Here  ^?e,p.  can  be  ot,tained  from  the  gradients  of  the  system 
nodal  point  displacements  vector  W.  These  in  turn  can  be  calcu¬ 
lated  from  the  equation  of  equilibrium  (Eq.  5.29): 


K(p)W  =  L 


K  W  +  K  ( p )  W  =  L  =  0 


(7.20) 

(7.21) 


K  ( p )  W 

«  ~  ~,P. 


-  K  W 
«  »  P  j  ~ 


(7.22) 


W  can  be  solved  from  equation  (7.20)  and  W  can  be  solved 

~'Pj 

from  equation  (7.22).  Using  the  rectangular  matrix  Xg  defined 
in  section  5.2,  we  get; 


' 

w 

fT  1 

xx,  Pj 

ixx 

W 

fT 

yy/pj 

* 

~yy 

w 

fT 

xyf  p  j 

(x,y) 

~xy 

T  X  W 
~e~e~ , p . 


(7.23) 


(x,y) 


The  calculation  of  K  was  discussed  in  connection 

~'pj 

with  the  flutter  constraint  gradients. 
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7.4  Gradients  of  the  Thickness  Constraints 


From  equation  (5.32)  we  get 


Cti,pj  ”  6ij 


where  is  the  Kronecker's  delta. 


(7.24) 


CHAPTER  8 


SEMI- INFINITE  PANEL  OPTIMIZATION 

The  problem  of  semi-infinite  panel  optimization  with 
flutter  constraint  is  considered  first  since  there  are  control 
theory  and  finite  element  solutions  to  this  problem  available 
for  comparison  in  literature  (references  7  and  8) .  First,  the 
computer  routines  to  determine  acr  are  developed,  then  the 
routines  to  calculate  the  flutter  gradients  based  on  equations 
of  section  7.2  are  coded.  The  results  are  checked  against 
gradients  calculated  numerically .  Once  these  routines  are 
established,  a  gradient  projection  routine  is  applied  to  opti¬ 
mize  the  panel  with  only  a  flutter  constraint.  This  routine 
had  already  been  developed  and  applied  to  problems  involving 
simple  algebraic  merit  functions  of  a  few  variables  with  non¬ 
linear  constraints  whose  gradients  were  easy  to  evaluate. 
Different  finite  element  combinations  are  tried  in  order  to 
find  the  best  representation  from  the  point  of  view  of  optimi¬ 
zation.  Later,  the  multiple-constraint  gradient  projection 
routine,  that  was  developed  primarily  for  the  delta  wing  pro¬ 
blem  is  applied  to  the  best  finite  element  representation  to 
study  the  effect  of  damping  on  optimization. 

8.1  Structural  Model 

The  following  assumptions  define  the  structural  model: 

1)  The  panel  has  length  t,  is  simply  supported  at 
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both  ends  in  the  direction  of  the  air  flow  and  is  free  of  in¬ 
plane  forces. 

2)  The  panel  is  of  infinite  dimension  in  the  direc¬ 
tion  perpendicular  to  the  air  flow,  i.e.  the  problem  is  one 
dimensional. 

3)  The  panel  is  of  sandwich  construction. 

4)  The  effect  of  any  air  below  the  panel  can  be 

neglected. 

Figure  8.1  shows  the  cross  section  of  the  structural 

model. 


Fig.  8.1  Semi-Infinite  Simply  Supported  Sandwich 
Panel  in  Supersonic  Flow 
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8.2  Analysis  to  Determine  a 
_ _ _ cr 

Five  constant  thickness  beam  elements  are  used  across 
the  span  l  to  represent  the  panel.  Mass,  stiffness,  aerody¬ 
namic  and  damping  matrices  for  a  constant  thickness  sandwich 
beam  element  are  given  in  appendix  A.  Span  length  l  is  used  to 
scale  the  finite  element  dimensions.  a  is  determined  for 

C  1C 

different  damping  values,  9ao* 

In  Table  8.1  first  two  eigenvalue  pairs  of  the  C 

~  o 

matrix  are  presented  for  increasing  a  values  and  for  different 

o 

damping  parameters,  g  .  The  starting  a  values  are  estimated 

ao  o 

from  a  graph  developed  by  Houbolt  (Ref.  23)  which  is  reproduced 
.  2 

in  reference  3.  There  is  a  ir  factor  difference  between  the 

damping  parameter  g^  used  by  Houbolt  and  the  damping  parameter 

g  used  in  this  study  such  that: 

^ao  J 

^ a  -  ^ao/n^  (8.1) 

Numbers  in  Table  8.1  exhibit  the  flutter  phenomenon 

clearly.  It  can  be  seen  that  for  small  damping  values  the  aQ 

for  which  the  imaginary  parts  become  equal  is  very  close  to  a  , 

CJT 

whereas  for  high  damping  values  otcr  is  much  higher  than  the  aQ 
for  which  the  frequencies  merge  (as  illustrated  in  figure  5.1). 
In  Table  8.2  otcr  values  and  flutter  frequencies  are  given  for 
all  of  the  damping  cases  considered.  In  Figure  8.2,  the  acr 
values  versus  <3aQ//v2  are  plotted  against  the  curve  given  by 
Houbolt,  (His  solutions  are  obtained  by  solving  the  differen¬ 
tial  equation  of  flutter) .  This  shows  that  the  a  values 

C  XT 
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TABLE  8.1  CRITICAL  EIGENVALUES  FOR  INCREASING  a  VALUES 

o 


ao 

h  =  <Bi  ?i  “i’/uo 

x2  =  (B±  Ti  “2)/wo 

£ 

0 

II 

o 

• 

o 

341.25 

342.50 

*  342.901 

343.75 

0.000  +i  31.20 

0.000  +i  31.81 

0.000  +i  32.38 

+  0.804  +i  32.40 

0.000  +i  33.45 

0.000  +i  32.92 

0.000  +i  32.38 

-  0.804  +i  32.40 

g  _  =  2.0 
^ao 

341.25 

343.75 

*  344.214 

345.00 

-  1.000  +i  31.19 

-  0.196  +i  32.39 

0.000  +i  32.40 

+  0.264  +i  32.43 

-  1.000  +i  33.43 

-  1.804  +i  32.39 

-  2.000  +i  32.40 

-  2.264  +i  32.43 

g  =  1  .Off2  g_  /0  =  4.935 

^ao  ^ao/2 

361.25 

367.50 

*  374.570 

400.00 

-  1.168  +i  32.59 

-  0.579  +i  32.80 

-  0.001  +i  33.03 

+  1.651  +i  33.87 

-  8.702  +i  32.59 

-  9.290  +i  32.80 

-  9.869  +i  33.03 

-11.52  +i  33.87 

gao  =  2.0»2  gao/2  =  9.87 

457.50 

463.75 

*  468.594 

477.50 

-  0.414  +i  34.83 

-  0.178  +i  35.06 

0.000  +i  35.23 

+  0.315  +i  35.54 

-19.33  +i  34.83 

-19.56  +i  35.06 

-19.74  +i  35.23 

-20.05  +i  35.54 

cr 


*Indicates  a 


TABLE  8.2  PANEL  FLUTTER  BOUNDARIES  AND  FLUTTER 
FREQUENCIES  FOR  DIFFERENT  DAMPING  VALUES 


g 

^ao 

a 

cr 

X  ~  =  +  i  (d),  ,  ) 

f  f/wo 

0.0 

342.901 

+i  32.38 

1.0 

343.230 

+i  32.38 

N> 

• 

O 

344.214 

+i  32.40 

1 .  Oir2 

374.570 

+  i  33,03 

1.5t r2 

414.375 

+i  33.95 

2.0tt2 

468.594 

+i  35.54 

Fig.  8.2  Flutter  Boundaries  for  Different  Damping  Values 


75 


obtained  by  using  five  constant  thickness  beam  elements  are 
slightly  lower  but  in  good  agreement  with  the  analytical  solu¬ 
tions  . 

8.3  Optimization  Results 

The  first  set  of  optimization  results  were  obtained 
by  using  a  gradient  projection  subroutine  which  can  only  handle 
one  constraint.  It  is  based  on  equations  (6.39  -  43).  For 
this  set  of  results  only  the  flutter  constraint  was  considered 
and  it  was  kept  active  all  the  times.  Various  methods  of  R 

q 

selection  and  various  convergence  criteria  were  used  with 
different  success  in  terms  of  number  of  iterations  and  the 
accuracy  of  the  final  results.  However  the  main  purpose  of 
this  initial  study  was  to  determine  a  proper  finite  element 
model  to  get  the  optimum  shapes  comparable  to  other  methods 
(mainly  control  theory  method  used  by  Weisshaar,  in  ref.  7)  . 

The  first  optimization  attempt  was  made  using  the 
same  model  used  for  calculating  the  acr  values,  i.e.,  five  con¬ 
stant  thickness  beam  elements  across  the  span  with  a  mass  ratio 
associated  with  each  element.  Using  equation  (A. 2) 

the  merit  function  for  this  model  can  be  written  as 


F(p)  =  l  P.j 

i=l 


(8.2) 


with  the  gradient 


VF  (p) 1  =  {1,  1,  1,  1,  1} 


(8.3) 


ao 


=  l.Oir  was  used  to- 


A  damping  parameter  of  g 
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gether  with  aQ  =  375.00  which  is  slightly  (0.11%)  higher  than 

the  a  =  374.57  (Table  8.1)  for  this  damping.  n  =  0.80  was 
cr 

used  for  the  ratio  of  initial  cover  skin  mass  to  total  mass. 

The  initial  flutter  gradients,  VC^(po),  were  checked 

against  numerical  gradients  which  were  obtained  by  calculating 

-5 

the  flutter  eigenvalue  with  a  small  increment  (e  =  10  )  to 

each  design  variable  successively.  The  agreement  was  found 

to  be  extremely  good.  During  the  optimization,  the  gradients 

were  calculated  for  all  p.  and  the  symmetry  property  was  used 

3Cf  3Cf 

to  check  the  accuracy  of  the  gradients  (i.e.  —  =  -r— —  , 

3C.  3Cf  P1  _  P5 

— 2.  =  — .  See  ref.  4  for  the  proof  that  the  optimum  shape 

0Q2  °P4 

is  symmetrical  about  mid  span) . 

The  initial  and  final  values  for  the  merit  function, 
the  design  parameters  (p^  and  p^  are  omitted  because  of 
symmetry) ,  and  the  flutter  eigenvalue  pair  are  given  in  table 
8.3.  The  results  are  also  plotted  in  figure  8.3.  It  can  be 
seen  that  thickness  constraints  never  come  into  the  picture. 

The  other  point  of  interest  is  that  the  flutter  frequency  of 
the  optimum  shape  is  1.62%  higher  than  the  initial  value.  This 
model  produced  4.96%  structural  weight  reduction  which  means 
3.97%  reduction  in  total  weight.  These  results  are  very  close 
to  the  results  presented  by  Craig  (Ref.  8)  using  5-constant 
thickness  elements  and  zero  damping. 

The  crudeness  of  the  optimum  shape  for  constant 
thickness  element  suggested  that  linearly  tapered  elements  be 
used.  Therefore  the  same  problem  was  solved  using  5  tapered 
elements  across  the  span.  The  element  mass,  stiffness,  aero- 
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TABLE  8.3  PANEL  OPTIMIZATION  -  5  CONSTANT  THICKNESS  ELEMENTS 
(NO  THICKNESS  CONSTRAINT) 


n  =  0.80 


1.000  1.000  1.000 


g  =  l.Oir" 
yao 


F(p) 


INITIAL  5.000 


FINAL 


TABLE  8.4  PANEL  OPTIMIZATION  -  5  TAPERED  ELEMENTS 
(NO  THICKNESS  CONSTRAINT) 


ao 

~ 

375.00 

0.0319 

+i 

33.041 

0.0319 

+i 

33.576 

g  =  l.Oir 
^ao 


n  =  0.80 


a  =  375.00 
o 


F  (p) 


INITIAL  5.000  1.000  1.000  1.000  0.0319  +i  33.041 


FINAL  4.546  0.0795  1.2457  0.9479  0.0319  +i  33.268 

_ _ . _ _ _ _ _ _ _  _  -  ■  - -  ■  - - 


TABLE  8.5  PANEL  OPTIMIZATION  -  6  TAPERED  ELEMENTS 
(NO  THICKNESS  CONSTRAINT) 
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dynamic  and  damping  matrices  for  tapered  sandwich  beam  element 
are  in  appendix  B.  This  time  a  mass  ratio  was  associated 
with  each  nodal  point  which  increased  the  number  of  design 
variables  from  5  to  6.  The  merit  function  for  this  model  is 
(using  equation  B.2); 


F(p) 


1 

2P1 


+  l  pi 

i=2 


(8.4) 


and  the  merit  function  gradient  is; 

VF(p)T  =  {|,  1,  1,  1,  1,  j>  (8.5) 


The  results  are  given  in  table  8.4  and  figure  8.4. 
This  time  the  increase  in  flutter  frequency  is  0.69%.  The 
structural  weight  reduction  is  9.08%  and  the  total  weight 
reduction  is  7.27%,  which  show  a  considerable  improvement  over 
constant  thickness  element  representation.  However  when  the 
optimum  shape  in  figure  8.4  is  compared  with  the  optimum  shapes 
obtained  by  control  theory  methods  (Ref.  7,  pp.  129,  132,  133) 
the  inability  of  5  tapered  element  model  to  reproduce  the  dip 
at  the  center  became  obvious.  This  suggested  the  use  of  6 
tapered  elements. 

2 

The  same  problem  (i.e.  gaQ  =  l.Oir  ,  aQ  =  375.00, 
ri  =  0.80)  was  also  solved  by  using  six  tapered  elements  with 
merit  function; 

1  6  1 

F(p)  =  +  I  P±  +  jP?  (8.6) 

2.  2 


and  merit  function  gradient 


i 


VF(p) 


{1/2,  1,  1,  1,  1,  1,  1/2} 


(8.7) 
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The  results  are  in  table  8.5  and  figure  8.5.  It  can 
be  seen  that  this  model  gave  the  required  results  as  far  as  the 
optimum  shape  is  concerned.  However,  it  also  indicated  the 
necessity  of  using  some  sort  of  thickness  constraints.  Without 
increasing  the  number  of  constraints  two  different  approaches 
were  tried  with  some  success: 

1)  Once  the  condition  p.  <  p  .  is  encountered  p. 

i  mm  l 

was  replaced  by  Pm^n  and  this  design  variable  was  excluded  from 
the  rest  of  the  optimization  process. 

2)  Every  time  the  condition  p.  <  p  .  is  encountered 

l  min 

was  replaced  by  Pm^n  but  kept  as  a  design  variable. 

Both  of  these  methods  were  effective  in  applying  the 
thickness  constraints.  However  both  had  undesirable  aspects: 
Method  1)  reduced  the  size  of  the  optimization  problem,  which 
meant  solving  not  the  original  problem  but  a  projection  of  it  on 
a  subspace.  The  optimum  design  obtained  this  way  was  not  the 
true  optimum.  Method  2)  interfered  with  the  calculations  of 
the  gradient  projection  method  which  either  meant  more  iter- 
actions  or  no  convergence. 

After  the  gradient  projection  routine,  which  could 
handle  multiple  constraints,  was  developed  for  the  delta  wing 
problem,  it  was  applied  to  the  panel  problem.  The  same  struc¬ 
tural  model  was  used  (i.e.  six  tapered  elements),  Pm^n  =  0.10 
and  n  »  0.70  were  chosen  to  form  a  better  comparison  with  the 
optimum  shape  given  by  Weisshaar  (Ref.  7,  p.  129) .  Each  thick¬ 
ness  constraint  was  handled  as  a  separate  constraint  according 


375.00 
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to  the  equations  (5.32)  and  (7.24).  Different  damping  values 

were  used  to  investigate  the  effect  of  damping  on  the  optimized 

2 

shape  and  weight  reduction.  gaQ  =  O.OItt  was  used  instead  of 
of  zero  damping.  This  was  done  to  avoid  trouble  in  calculating 
the  eigenvectors,  since  at  flutter  there  are  repeated  eigen¬ 
values  for  zero  damping  and  repeated  eigenvalues  cause  problems 
with  inverse  iteration. 

The  results  of  importance  are  given  in  tables  8.6-9. 

2 

In  figure  8.6  the  optimum  shape  for  gaQ  =  O.OItt  is  plotted 

against  the  control  theory  result  given  by  Weisshaar  who  uses 

zero  damping  and  the  same  Pm^n  and  n .  This  plot  shows  that  the 

six  tapered  element  model  is  very  efficient  in  approximating  the 

optimum  shape.  In  figure  8.7  all  the  optimum  shapes  for 

different  dampings  are  plotted. 

The  results  show  that  the  damping  practically  has  no 

2  2 

effect  on  the  optimum  shape  up  to  g  =  I.Ott  .  For  g  =  1.5tt 

060  UU 

the  optimum  shape  is  somewhat  changed  but  still  similar.  How- 

2 

ever  for  g  =  2.0tt  the  optimum  shape  is  drastically  effected, 
ao 

The  results  also  show  that  the  effect  of  damping  on  weight 

2 

reduction  is  favorable,  however  slight,  up  to  gaQ  =  1.5ir  .  The 
structural  and  total  weight  reductions  for  different  models  and 
for  different  damping  values  are  given  in  table  8.10.  Percen¬ 
tages  taken  from  Weisshaar' s  work  (Ref.  7)  belong  to  the  optimum 
shape  reproduced  in  figure  8.6. 

The  increases  in  the  flutter  frequency  from  initial  to 

2 

final  design  are  4.45%,  6.32%,  8.72%  and  63.8%  for  9ao?  O.Olir  , 

2  2  2 
I.Ott  ,  1.5tt  ,  and  2.0-rr  respectively. 


TABLES  8.6-9  PANEL  OPTIMIZATION  WITH  DIFFERENT 
DAMPING  VALUES  -  6  TAPERED  ELEMENTS 


(THICKNESS  CONSTRAINTS) 


gaQ  =  0.01*  n  =  0.70 


F(p) 


p  .  =  0.1  a  =  343.1375 

Kmin  o 


INITIAL  6.000  1.000  1.000  1.000  1.000  0.0461  +i  32.404 


FINAL  5.146  O.lOOa  1.0853  1.3878  0.1000  0.0477  +i  33.845 


g  =  l.O*^ 

n  =  0.70 

p  .  =0.1 

a  =  375.000 

^ao 

o 

F  (p) 


Xf 

0.0021 

+i 

33.068 

0.0024 

+i 

35.161 

g  =  1.5* 
^ao 


n  =  0.70 


p  .  =  0.1  a  =  414.375 

Kmin  o 


F(p) 


INITIAL  6.000  1.000  1.000  1.000  1.000  0.0016  +i  33.954 


FINAL  4.852  0.1000  1.1287  1.1971  0.1000  0.0017  +i  36.917 


g  =  2.0* 
^ao 


n  =  0.70 


p  .  =  0.1  a  =  469.625 

min  o 


F  (p) 


INITIAL  6.000  1.000  1.000  1.000  1.000  0.0016  +i  35.310 


FINAL  2.551  0.1000  0.8490  0.3263  0.1000  0.0026  +  57.830 
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2 

Damping  values  over  gaQ  =  I.Ott  are  rare  for  most  lift¬ 
ing  surfaces?  however,  they  are  possible,  especially  for  members 
where  a  special  viscous  core  material  is  used  for  noise 
suppression. 

To  give  an  idea  of  the  optimization  process,  step  by 

2 

step  results  are  given  in  table  8.11  for  the  gaQ  =  I.Ott  case. 

The  step  size,  R^,  pattern  for  this  run  was  predetermined,  using 
earlier  experience  with  single  constraint  gradient  projection 
algorithm.  In  figure  8.8  merit  function  and  flutter  constraint 
values  are  plotted  against  step  sizes  (positive  constraint  values 
indicate  constraint  violation) .  The  option  of  reducing  step  size 
when  a  new  thickness  constraint  is  encountered  (Equation  6.38) 


has  not  been  used  for  these  runs . 


TABLE  8.11  STEP  BY  STEP  OPTIMIZATION  -  6  TAPERED  ELEMENTS 


Fig.  8.8  Progress  of  Panel  Optimization  -  6  Tapered  Elements 
(gao  =  I.Ott2,  n  =  0.70,  pm.n  =  0.10,  aQ  =  375.00) 


CHAPTER  9 


CANTILEVERED  DELTA-WING  OPTIMIZATION 

The  results  of  semi-infinite  panel  optimization  showed 
that  the  aeroelastic  optimization  approach  presented  in  chapters 
2,  3,  4,  5,  6  and  7  is  basically  feasible.  The  second  problem 
attempted  here  is  the  minimum-weight  design  of  a  cantilevered 
delta  wing  with  flutter,  stress  and  thickness  constraints.  This 
problem  is  two  dimensional  and  considerably  larger  in  size  com¬ 
pared  to  the  panel  problem.  The  same  structural  model  is  used 
in  the  optimization  problem  with  three  different  combinations 
of  constraints:  1)  Flutter  and  thickness  constraints  2)  Stress 
and  thickness  constraints  3)  Flutter,  stress  and  thickness  con¬ 
straints.  The  multiple-constraint  gradient  projection  routine 
which  was  discussed  in  Chapter  6  was  first  applied  to  the  opti¬ 
mization  problems  including  the  flutter  constraint. 

9.1  Structural  Model 

The  structural  model  is  a  cantilevered  delta  wing  with 
aspect  ratio  1.  It  is  assumed  to  be  of  sandwich  construction 
with  constant  core  thickness  d  and  with  variable  skin  thickness 
t.  The  air  flow  is  assumed  to  be  parallel  to  the  cantilevered 
edge  and  at  both  sides  of  the  wing.  Figure  9.1  shows  the  plan 
view  of  the  delta  wing.  The  finite  element  lay-out,  numbering 
of  the  nodal  points  and  the  elements  are  also  shown  in  figure 
9.1.  The  length  of  the  clamped  edge  is  L,  which  is  used  to 
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scale  all  finite  element  dimensions. 

A  high-precision  triangular  plate  bending  element  deve¬ 
loped  by  Cowper  (Ref.  22)  and  others  is  used  for  the  finite  ele¬ 
ments.  A  design  variable  is  associated  with  each  nodal  point. 
The  cover  skin  thickness  variation  within  an  element  is  assumed 
to  be  linear.  The  equations  for  the  generation  of  element 
stiffness  and  mass  matrices  for  a  constant  thickness  plate  ele¬ 
ment  are  given  in  reference  22.  These  equations  are  modified  to 
obtain  the  equations  for  a  sandwich  plate  with  linear  cover  skin 
thickness  distribution.  The  equations  for  the  generation  of 
element  aerodynamic  and  damping  matrices  are  derived.  The  ele¬ 
ment  characteristics,  and  the  equations  for  the  generation  of 
element  mass,  stiffness,  aerodynamic  and  damping  matrices  are 
given  in  appendix  C.  The  actual  element  matrices  are  generated 
numerically  inside  the  computer. 

There  are  six  nodal  point  degrees  of  freedom,  namely 
the  transverse  deflection  and  its  first  and  second  derivatives 

(W,  W  ,  W  ,  W  ,  W  ,  W  ) .  The  system,  therefore,  has  60  de- 
'  '  x'  yf  xx'  xy  y  J 

grees  of  freedom  before  the  boundary  conditions  are  applied. 

Prescribed  boundary  conditions  are  applied  for  the  nodal  points 

on  the  clamped  edge,  which  are  1,  2,  4  and  7  (Figure  9.1).  The 

fixed-edge  condition  implies  that  W  =  Wy  =  0  along  the  x-axis, 

from  which  we  also  get  W  =  W  =  W  =  0.  Therefore  nodal 

x  xy  xx 

points  1,  2,  4,  and  7  will  only  have  Wyy  freedom.  Enforcement 
of  the  prescribed  boundary  conditions  reduces  the  system  degrees 
of  freedom  to  40.  These  boundary  conditions  are  essential  for 
any  type  of  structural  analysis,  and  they  can  be  assumed  suffi- 
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cient  for  a  load-deflection  or  a  dynamic  (natural  frequencies, 
flutter  boundaries  etc.)  type  analysis.  However  the  second 
derivatives  of  the  transverse  deflection  (Wxx,  WXy,  W  )  ob¬ 
tained  this  way  for  nodal  points  on  free  edges  (nodal  points  3, 

6,  10,  9  and  8  at  figure  9.1)  will  not  satisfy  the  free  edge 
natural  boundary  conditions  for  plate-bending  moments  (bending 
moments  and  second  derivatives  are  related  through  equations 
3.2).  For  the  structural  model  used  in  optimization,  in  order 
to  obtain  more  realistic  stress  constraints  for  free  edge  nodal 
points,  these  degrees  of  freedoms  are  constrained  to  satisfy 
the  free  edge  conditions.  These  natural  boundary  conditions 
reduce  the  number  of  system  degrees  of  freedom  to  29  (the  alge¬ 
braic  details  for  the  application  of  boundary  conditions  are  pre¬ 
sented  in  appendix  C) . 

The  40  degree  of  freedom  system  is  used  to  obtain  the 
natural  frequencies  and  flutter  boundaries  of  the  delta  wing 
which  are  compared  to  the  results  given  in  literature  (Refs.  22, 
6) .  The  29  degree  of  freedom  system  is  also  used  in  the  same 
calculations  to  see  the  effect  of  additional  constraining  on 
natural  frequencies  and  flutter  boundaries.  Only  the  29  degree 
of  freedom  system  is  used  for  the  optimization  procedures. 

Poisson's  ratio  v  =  0.3  and  skin  mass  ratio  n  =  0.80  is 
used  in  all  the  calculations  presented  below. 

9.2  Natural  Frequencies 

Natural  frequencies  are  calculated  to  check  the 
accuracy  of  the  derivation  of  system  mass  and  stiffness  matrices. 
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TABLE  9.1  NATURAL  FREQUENCIES  OF  THE  DELTA  WING 


Mode  No. 

Natural  Fre¬ 
quencies  Given 
by  Cowper  et  al 
(Ref.  22) 

Natural  Fre¬ 
quencies  of 
40  D.O.F. 
System 

Natural  Fre¬ 
quencies  of 
29  D.O.F. 
System 

% 

Difference 

■■ 

36.6419 

36.6419 

37.0083 

1.0% 

2 

139.3265 

139.3265 

142.7298 

2.4% 

3 

194.1408 

194.1408 

198.3996 

2.2% 

4 

333.829 

333.829 

338.277 

1.3% 

5 

455.374 

455.374 

459.011 

0.8% 

6 

593.238 

593.238 

617.275 

1 

4.0% 

7 

671.342 

1 

671.342 

727.063 

8.4% 

8 

811.612 

811.612 

848.693 

4.5% 

9 

969.650 

969.650 

1030.104 

6.2% 

10 

1126.614 

1126.614 

1171.494 

■SB 
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The  equation; 

M”1  K  W  =  W  (9.1) 

is  used  for  this  purpose.  In  reference  22  the  first  10  natural 
frequencies  are  given  for  a  cantilevered  triangular  plate  with 
aspect  ratio  1  and  Poisson's  ratio  v  =  O.3..  The  finite  element 
lay  out  is  exactly  the  same.  There,  actual  material  and  dimen¬ 
sional  properties  are  used.  Since  the  matrices  derived  here 
are  non-dimensional,  a  comparison  is  obtained  by  scaling  the 
natural  frequencies  obtained  here  with  a  factor  such  that  the 
first  natural  frequencies  are  identical.  Both  40  and  29  degree 
of  freedom  systems  are  used  in  these  calculations.  The  scaling 
factor  is  obtained  with  respect  to  the  first  natural  frequency 
of  the  40  degree  of  freedom  system.  The  results  are  presented 
in  Table  9.1.  The  natural  frequencies  obtained  by  the  40  de¬ 
gree  of  freedom  system  are  exactly  the  same  as  the  ones  given 
by  the  authors  of  reference  22.  The  natural  frequencies  ob¬ 
tained  by  29  degree  of  freedom  system  are  slightly  higher  be¬ 
cause  of  the  extra  constraining.  On  the  average,  the  first 
five  frequencies  are  1.5%  higher,  and  the  last  five  frequencies 
are  5.4%  higher. 

9 . 3  Flutter  Boundaries 

In  reference  6  Olson  gives  the  flutter  boundaries  for 
cantilevered  delta  wings,  for  different  sweepback  angles  A 
(Figure  9.1)  and  for  different  finite  element  grids.  He  uses 
the  same  aerodynamics  with  zero  damping  and  the  same  finite  ele¬ 
ments.  For  A  =  45°  and  for  3x3  grid  the  flutter  boundary  is 


TABLE  9.2  CRITICAL  EIGENVALUES  FOR  40  D.O.F.  SYSTEM 


ll  "  (01  +x  “lVwo 


2  =  (B2  +i  “2}/«o 


-  0.1234  +i  6.172  -  0.1234  +i  23.473 


70.0 

-  0.1234  +i  14.196 

-  0.1234  +i  21.418 

86.245 

-  0.1213  +i  18.567 

-  0.1255  +i  18.573 

86.3 

+  0.0688  +i  18.572 

-  0.3156  +i  18.572 

86.4 

+  0.2583  +i  18.577 

-  0.5051  +i  18.577 

88.0 

+  1.1949  +i  18.653 

-  1.4416  +i  18.653 

TABLE  9.3  CRITICAL  EIGENVALUES  FOR  29  D.O.F.  SYSTEM 

(g  =  0.2468) 


86.35 


91.00 


95.00 


95.50 


*95.60 


107.00 


‘1  *  (ei  “l>/MO 


-  0.1234  +i  16.275 


-  0.1234  +i  17.281 


-  0.1234  +i  18.682 


-  0.1234  +i  19.159 


+  0.0754  +i  19.352 


2  =  (e2  +i  “2)/«o 


-  0.1234  +i  21.535 


-  0.1234  +i  20.976 


-  0.1234  +i  19.963 


-  0.1234  +i  19.536 


-  0.3221  +i  19.352 


+  2.8007  +i  19.914  -  3.0475  +i  19.914 


♦Indicates  The  Value  Used  For  Optimization 
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given  as  a  =  86.254.  Here  the  eigenvalues  of  the  C  matrix 

cr  ~o 

are  calculated  for  increasing  aQ  values.  A  damping  value  of 

g  „  =  0.025tt^  (=  0.2468)  is  used.  In  table  9.2  the  first  two 

ao 

eigenvalue  pairs  for  the  40  degree  of  freedom  system  are  pre¬ 
sented.  a  is  determined  roughly  to  be  a  =86.3.  The  slight 
difference  can  be  explained  with  the  existence  of  small  damping 

which  does  tend  to  increase  the  a  .  Table  9.3  is  for  the  29 

cr 

degree  of  freedom  system  where  a  is  determined  roughly  to  be 

C1T 

a  =95.6.  This  is  a  10.8%  increase  from  the  a  for  the  40 
cr  cr 

degree  of  freedom  system.  The  calculation  of  the  flutter  boun¬ 
dary  is  a  stability  problem  and  the  elimination  of  any  existing 
free  edge  moments  could  explain  the  rather  considerable  improve¬ 
ment  on  the  flutter  boundary  (for  the  40  degree  of  freedom  sys¬ 
tem  the  free  edge  moments  are  not  zero) .  The  increase  in  the 
flutter  frequency  is  4.2%. 

Using  a  proper  scaling  factor  the  flutter  frequency  for 
the  40  degree  of  freedom  system  is  found  to  be  in  very  good 
agreement  with  the  results  given  by  Olson  (Ref.  6) . 

9.4  Optimization  Results 

Using  equation  (C.64)  given  in  appendix  C  for  the  weight 
of  a  triangular  element,  the  merit  function  can  be  derived  to 
be? 

F  ( p )  =  +  3p2  t  3p^  -^4  6pij  3Pg  P y  3pg 


+  3p9  +  p1Q 


(9.2) 
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From  which  the  merit  function  gradients  are  obtained  as; 

VF(p)T  =  {1,  3,  3,  3,  6,  3,  1,  3,  3,  1}  (9.3) 

These  expressions  are  used  in  obtaining  all  of  the 
following  optimization  results. 

9.4.1  Optimization  with  Flutter  rnd  Thickness  Constraints 

The  first  results  presented  here  are  for  the  optimi¬ 
zation  of  the  delta  wing  with  flutter  and  thickness  constraints. 
The  gradients  of  the  flutter  constraint  are  checked  numerically 
(as  described  for  the  panel  problem)  and  found  to  be  extremely 
accurate.  The  29  degree  of  freedom  system  is  used  in  obtaining 
the  system  matrices.  Other  parameters  for  the  optimization  are: 

g  =  0.025tt^,  a  =95.6  (From  table  9.3)  and  n  =  0.80.  The 

ao  cr 

small  damping  value  is  used  to  separate  critical  eigenvalues 
at  flutter  which  increases  the  efficiency  in  calculating  the 
corresponding  eigenvectors . 

In  table  9.4  the  step  sizes,  the  design  variables,  the 
merit  function,  the  flutter  eigenvalue,  the  flutter  constraint 
and  the  active  thickness  constraints  are  presented  for  21  de¬ 
sign  steps.  The  flutter  constraint  is  kept  active  for  all  de¬ 
sign  steps  and  it  is  defined  as  =  (8-  -  ftf0)/w0  since  6f0 

corresponding  to  a  =95.6  was  slightly  greater  than  zero. 

These  results  were  obtained  in  four  different  computer  runs; 
the  first  three  runs  for  six  design  steps  each  and  the  last  run 
for  three  design  steps.  On  the  University  of  Texas  at  Austin 
CDC  6600  each  design  step  took  approximately  21  seconds.  The 
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TABLE  9.4  DELTA  WING  OPTIMIZATION  WITH  FLUTTER 
AND  THICKNESS  CONSTRAINTS 


Given  R_ 


• 

o 

0.4 

0.4 

0.3 

0.3 

0.4 

- 1 - 

0.4 

- , - 

0.4 

- 1 - 

0.290 

*"  '"1  *  ■  ■■ 

0.3 

The. 

Con. 


1.000 

0.9684 

0.9414 

0.9222 

0.9169 

0.8977 

1.000 

0.8933 

0.7977 

0.7214 

0.6881 

0.6218 

1.000 

0.9476 

0.9202 

0.9225 

0.9376 

0.8769 

1.000 

0.8604 

0.7233 

0.5924 

0.5086 

0.4180 

1.000 

- 1 

0.8473 

0.6863 

0.5072 

0.3390 

0.1118 

1.000 

1.0985 

1.1789 

— - 1 

1.2272 

1.2037 

1.0982 

1.000 

0.9748 

0.9506 

0.9292 

0.9162 

0.8901 

1.000 

0.8340 

0.6664 

0.5015 

.  | 

0.3979 

0.3067 

1.000 

0.7697 

0.5287 

0.2773 

0.1000 

0.1 

1.000 

0.8947 

0.7876 

0.6809 

0.6129 

0.5730 

F(p) 

27.000 

24.132 

21.243 

18.302 

15.988 

13.296 

^f/uo 

0.0754 

-0.0764 

-0.0778 

-0.0793 

-0.0818 

-0.0903 

wf/wo 

19.352 

18.618 

17.827 

17 .035 

16.624 

16.393 

cf 

0.00 

-0.1518 

-0.1532 

-0.1547 

-0.1572 

-0 .1657 

-7.17 


-7.22 


-7.35 


-8.97  -7.28 


gao  0 . 025tt  ,  acr  -  95.6,  n  -  0.80,  Cf  -  (8fg  -  3fo)/u 
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TABLE  9.4  (Continued)  DELTA  WING  OPTIMIZATION  WITH  FLUTTER 

AND  THICKNESS  CONSTRAINTS 


0.3 

warn 

0.2 

CM 

• 

O 

0.2 

0.2 

R 

q 

m 

0.2 

0  .164 

0.2 

_ 1 

0.139 

_ , _ , 

0.182 

6 

7 

8 

9 

10 

11 

a 

0.8796 

0.8470 

0.8189 

0.7749 

0.7401 

0.7125 

p2 

0.5721 

0.4775 

0.3977 

0.2800 

0.1973 

0.1000 

b 

0.8367 

0.7487 

0.6756 

0.5741 

0.5015 

0.4020 

a 

0.3804 

0.2941 

0.2257 

0.1422 

0.1000 

0.1 

B 

0.1000 

0.1 

0.1 

0.1 

0.1 

0.1 

B 

1.0860 

1.0159 

0.9601 

0.8956 

0.8456 

0.7408 

B 

0.8769 

0.8477 

0.8240 

0.7925 

0.7715 

0.7369 

P8 

0.2641 

0.1739 

0.1000 

0.1 

0.1 

0.1 

B 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

p10 

0.5594 

0.5295 

0.5032 

0.4611 

0.4247 

0.3990 

F(p> 

12.634 

11.254 

10.123 

8.904 

8.069 

7.077 

-0.0905 

-0.0906 

-0 .0908 

-0.0914 

-0.0916 

/  too 

16.330 

16.741 

17.120 

17.678 

18.171 

19.515 

cf 

-0.1659 

-0.1660 

-0.1622 

-0.1668 

-0.1670 

The. 

Con. 

5,9 

5,9 

5,8 

,9 

5,8 

,9 

4,5 

a\ 

00 

2,4 

8,9 

/  5 

AF/R 

-6.90 

00 

00 

• 

VO 

I 

-6.10 

o 

° 

j  • 

vo 

1 

-5.47 

q 

gao  =  0.025tt2,  acr  =  95.6,  n  =  0.80,  Cf  =  (6fq  -  Bfo>  A>0 
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TABLE  9.4  (Continued)  DELTA  WING  OPTIMIZATION  WITH  FLUTTER 

AND  THICKNESS  CONSTRAINTS 


gao  =  0 .025tt  ,  acr  =  95.6,  n  =  0.80,  Cf  =  (3f  -  3fo)/a)Q 
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TABLE  9.4  (Continued)  DELTA  WING  OPTIMIZATION  WITH  FLUTTER 

AND  THICKNESS  CONSTRAINTS 


AF/R„ 


Given  R 

q 

0.1 

0.05 

0.05 

0  .05 

R 

q 

0.1 

0.05 

0.05 

0.05 

Design 

No. (q) 

00 

19 

— 

20 

21 

Pi 

0.6006 

.  ■■■  ■ 

0.5988 

0.5907 

_ 

0.5766 

P2 

0.1 

0.1 

0.1 

p3 

0.2154 

0.2614 

0.2205 

0.2437 

P4 

0.1 

0.1 

0.1 

0.1 

P5 

0.1 

0.1 

0.1 

-  . 

p6 

0.5900 

0.6004 

0.5757 

- . -  -----  - 

0.5521 

p7 

0.5568 

0.5577 

0.5456 

0.5358 

P8 

0.1 

0.1 

mn 

0.1 

P9 

0.1 

0.1 

0.1 

0.1 

p10 

0.2551 

0.2388 

0.2354 

0.2021 

F(p) 

5.629 

5.781 

5.560 

5.502 

^f/wo 

2.897 

-0.0084 

2.311 

1.838 

wf/wo 

32.060 

24.514 

31.794 

27.754  | 

cf 

2.822 

-0.0838 

2.236 

1.762 

The. 

2,4,5 

2,4,5 

in 

nr 

CM 

_ - 

2,4,5 

Con. 

8,9 

8,9 

8,9 

""1  J 

8,9 

"  -  1  1 - 

-0.32 

+3.04 

o 

nr 

• 

nr 

1 

-1.17 

gao  =  °-025"2'  acr  =  95*6'  0  “  0*80,  Cf  =  (6fq  -  6f0)/w0 
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calculation  of  the  eigenvalues  and  the  left  and  right  hand 
eigenvectors  corresponding  to  the  flutter  eigenvalue  of  the  C 
matrix  takes  approximately  14  seconds. 

Table  9.4  reveals  that  the  first  thickness  constraint 
violation  occurs  at  the  4th  design  step  at  nodal  point  9  which 
is  the  reason  for  reduction  of  R3  (i.e.  norm  of  the  design 
change  vector  from  p3  to  p4)  from  0.30  to  0.290.  The  second 
thickness  constraint  is  encountered  at  the  6th  design  step  for 
nodal  point  5  which  causes  R5  to  be  reduced  from  0.3  to  0.091. 
Thickness  constraints  for  nodal  points  8,  4  and  2  are  encountered 
at  design  steps  8,  10  and  11  respectively.  Each  time  a  new 
thickness  constraint  is  encountered  there  is  a  reduction  in  the 
given  step  size  due  to  the  provision  discussed  in  Chapter  6 
(Equation  6.38). 

Until  the  11th  design  step  the  flutter  constraint  is 
closely  followed  within  the  feasible  region.  At  design  steps 
12,  14,  16,  18  and  20  the  flutter  occurs  involving  the  second 
and  third  modes.  To  clarify  this  phenomenon  the  first  three 
eigenvalue  pairs  are  given  in  table  9.5  for  several  optimiza¬ 
tion  steps.  It  can  be  seen  that  the  imaginary  parts  of  the 
eigenvalues  of  the  second  and  third  modes  come  closer  as  the 
optimization  progresses  until  the  12th  design  step  where  they 
become  close  enough  for  the  energy  transfer  necessary  for  the 
flutter  phenomenon. 

In  figure  9.2  the  merit  function  F(p)  and  flutter  con¬ 
straint  Cf  are  plotted  against  step  sizes.  It  is  interesting 
to  note  that  there  is  very  little  weight  reduction  after  11th  de- 
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TABLE  9.5  FIRST  THREE  EIGENVALUE  PAIRS 
DURING  OPTIMIZATION 


Design 
No.  (q) 

X1=(B 

“l'/ioo 

^2  ^2 

+i 

u2  ^ /wo 

A 

3=(B3  ?i  0,3) /M0 

0 

0.075 

+i 

19.35 

- 

0.322 

+i 

19.35 

- 

0.123 

+i 

38.03 

5 

- 

0.090 

+i 

16.39 

- 

0.192 

+i 

25.87 

- 

0.276 

+i 

42.30 

10 

- 

0.091 

+i 

18.17 

- 

0.209 

+i 

26.94 

- 

0.345 

+i 

39.29 

11 

- 

0.092 

+i 

19.51 

- 

0.208 

+i 

27.75 

0.387 

+i 

37.33 

12 

- 

0.099 

+i 

22.49 

+ 

2.527 

+i 

32.11 

- 

3.192 

+i 

32.49 

13 

- 

0.042 

+i 

24.03 

- 

0.214 

+i 

29.27 

- 

0.510 

+i 

34.74 

14 

- 

0.093 

+i 

25.71 

+ 

4.191 

+i 

31.75 

- 

4.914 

+i 

32.06 

15 

+ 

1.429 

+i 

26.76 

- 

1.746 

+i 

27.23 

- 

0.473 

+i 

35.36 

16 

- 

0.079 

+i 

25.71 

+ 

3.994 

+i 

31.82 

- 

4.741 

+i 

32.14 

17 

+ 

2.419 

+i 

27.11 

- 

2.754 

+i 

27.46 

- 

0.474 

+i 

35.46 

18 

- 

0.063 

+i 

24.54 

+ 

2.897 

+i 

32.06 

- 

3.655 

+i 

32.42 

19 

- 

0.008 

+i 

24.51 

- 

0.235 

+i 

29.58 

- 

0.556 

+i 

34.57 

20 

- 

0.030 

+i 

25.44 

+ 

2.311 

+i 

31.79 

- 

3.107 

+i 

32.22 

21 

+ 

1.838 

+i 

27.75 

- 

2.116 

+i 

28.27 

- 

0.550 

?i 

34.09 

Design  No. 
(Scaled  by  R  ) 

*  q 

9.2  Progress  of  Delta  Wing  Optimization  with  Flutter  and 
Thickness  Constraints 
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F(pq)  =  27.00  F(pf)  =  5.78 

SWR:  79%  TWR:  63% 


Fig.  9.3  Skin  Thickness  Contours  for  Minimum-Weight 
Design  with  Flutter  and  Thickness  Constraints 
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sign  which  is  in  the  feasible  region.  The  main  reason  to  carry 
the  optimization  beyond  this  point  was  to  see  the  behaviour  of 
the  flutter  constraint.  At  design  step  19  the  merit  function 
is  F(p,Q)  =  5.781  which  is  minimum  among  feasible  designs.  This 

~iy  i 

correspond  to  78.6%  structural  and  62.9%  total  weight  reduction. 
It  can  be  shown  that  the  weight  reduction  at  the  11th  design 
step  corresponds  to  94%  of  the  weight  reduction  achieved  at  the 
19th  design  step.  It  is  interesting  to  note  that  using  smaller 
step  sizes  and  greater  number  of  design  cycles  best  possible 
was  to  obtain  F(p)  =  5.700.  The  corresponding  thickness  dis¬ 
tribution  was  very  similar  to  the  design  vector  p-^  of  table  9.4. 

In  figure  9.3  the  plan  view  of  the  delta  wing  with  skin 
thickness  contours  corresponding  to  the  minimum-weight  design 
(P19)  is  presented. 

9.4.2  Optimization  with  Stress  and  Thickness  Constraints 

The  second  set  of  results  presented  here  are  for  the 

optimization  of  the  delta  wing  to  satisfy  stress  and  thickness 

constraints.  The  stress  constraints  are  placed  at  the  nodal 

points  where  the  constraints  can  be  evaluated  directly  from 

the  second  derivatives  contained  in  the  nodal  point  degrees  of 

freedom  (Equation  5.26).  It  was  assumed  that  the  delta  wing  is 

uniformly  loaded.  The  load  level  is  so  adjusted  that  the  nodal 

point  with  the  highest  combined  stress  (Equation  5.23)  is  at 

yield  for  the  initial  uniform  design.  The  stress  constraints 

are  scaled  so  that  C  =  -100  means  zero  stress  and  C  = 

s  s 

yield. 


0  means 
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The  significant  facts  about  22  design  steps  are  pre¬ 
sented  in  table  9.6.  These  results  were  obtained  in  a  single 
run  which  took  approximately  90  seconds  on  the  CDC  6600.  Each 
design  cycle  takes  little  less  than  4  seconds.  Successively 
smaller  values  of  e  (same  as  defined  in  Chapter  6)  were  used  to 

determine  critical  stress  constraints  (C  >  -e) .  The  values 

s 

employed  were:  e  =  25  for  the  first  6  design  steps,  e  =  15 
for  the  following  6  design  steps  and  e  =  5  after  the  12th  de¬ 
sign  step.  In  table  9.6  the  values  of  the  critical  stress  con- 
straints  are  given.  The  only  active  constraint  for  designs  0, 
1,  2,  and  3  is  the  stress  constraint  at  nodal  point  7.  The  re¬ 
duction  of  from  0.30  to  0.115  is  due  to  the  violation  of 
stress  constraint  at  nodal  point  8.  The  reductions  in  R.. ,  R„ , 
Rg  and  R-^q  are  due  to  the  violation  of  thickness  constraints 
at  nodal  points  9,  6,  2  and  3  respectively.  The  reduction  in 
the  number  of  active  stress  constraints  at  design  step  12  is 
due  to  the  decrease  in  the  value  of  e  from  15  to  5.  The  re¬ 
ductions  in  R-^2  /  and  R^  are  due  to  the  stress  constraint 

violations  at  nodal  points  8,  9  and  4.  The  reduction  in  R21  is 
due  to  the  thickness  constraint  violation  at  nodal  point  1. 

At  the  22nd  design  step  there  are  all  together  10  active  con¬ 
straints  (5  thickness,  5  stress  constraints)  and  that  is  equal 
to  the  number  of  design  variables.  This  indicates  that  any 
further  improvement  is  not  possible.  Although  the  merit  func¬ 
tion  is  slightly  lower  for  p19,  p22  is  a  better  design  with 
respect  to  the  value  of  the  stress  constraint  at  nodal  point  9. 
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TABLE  9.6  DELTA  WING  OPTIMIZATION  WITH  STRESS 
AND  THICKNESS  CONSTRAINTS 


Given  R, 


0.4  0.3 


0.4  ;  0.115 


H 

1.000 

0.9552 

0.9101 

0.8645 

0.8511 

0.7923 

m 

1.000 

0.8786 

0.7572 

0.6354 

0.6002 

0.4677 

1.000 

0.8767 

0.7545 

0.6337 

0.5996 

0.5074 

p4 

1.000 

0.9587 

0.9220 

0.8913 

0.8854 

0.8752 

- - 

m 

1.000 

0.7551 

0.5132 

0.2759 

0.2104 

0.2032 

n 

1.000 

0.8684 

0.7367 

0.6050 

0.5671 

0.4349 

I - 

H 

1.000 

1.0729 

1.1476 

1.2241 

1.2474 

1.3095 

1.000 

0.8477 

0.6909 

0.5286 

0.4793 

0.5662 

K 

1.000 

0.8638 

0.7262 

0.3783 

1.000 

0.9554 

0.9107 

0.8656 

0.8525 

0.7906 

f(p) 

27.000 

23.396 

19.810 

16.250 

15.242 

13.800 

Ill 


TABLE  9.6  (Continued)  DELTA  WING  OPTIMIZATION  WITH 
STRESS  AND  THICKNESS  CONSTRAINTS 


0.3  0.2 


0.3  0.187 


0.2 

CM 

« 

O 

CM 

• 

O 

0.079 

0.038 

0.2 

0.6952 


0.2519 


0.3576 


0.8553 


0.1946 


0.2325 


1.4066 


0.7088 


0.1000 


0.6882 


11.476 


9 


7,8 


8 

9 

0.6441 

0.6254 

0.1405 

0.1141 

0.2787 

0.2463 

0.8479 

0.8549 

0.1929 

0.1786 

0.1385 

0.1000 

1.4611 

1.4803 

0.7716 

0.8099 

0.1 

0.1 

0.6340 

0.6100 

10.728 


9 


4,7,8 


10.463 


6,9 


4,7,8 


3 


-8.70 


10 


0.6149 


0.1000 


0.2277 


0.8581 


0.1685 


0.1 


1.4907 


0.8288 


0.5963 


10.357 


2,6,9 


4,7,8 


-0.03 


-8.66 


-3.36  -2.81 


-2.36  -2.17 
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TABLE  9.6  (Continued)  DELTA  WING  OPTIMIZATION 
WITH  STRESS  AND  THICKNESS  CONSTRAINTS 


0.2 

0.1 

0.1 

0.1 

0.1 

0.1 

R 

q 

.003 

0.031 

0.1 

0.1 

0.099 

0.1 

12 

13 

14 

15 

16 

17 

B 

0.5343 

0.5290 

0.4726 

0.4160 

0.3600 

0.2976 

wm 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1000 

0.1 

0.1 

0.1 

0.1 

0.1 

SI 

0.8297 

0.8289 

0.7944 

0.7592 

0.7244 

0.6742 

E9 | 

0.1487 

0.1269 

0.1179 

0.1096 

0.10  22 

0.1175 

Hi 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

B 

1.5628 

- j 

1.5722 

1.6149 

1.6569 

1.6981 

1.7365 

p8 

0.8964 

0.8788 

-  •  -  - 

0.8978 

0.9159 

0.9330 

0.9226 

0.1 

0.1 

0.1 

0.1 

0.1 

0.1 

i 

0.5092 

0.5043 

0.4463 

0.3880  j  0.3300 

..  -J 

0.2880 

f(p) 

9.877 

9.690 

9.517 

— - 1 

9.344 

9.174 

9.017 

The. 

Con. 

2, 3,6, 9 

2, 3, 6, 9 

2, 3, 6, 9 

2, 3, 6, 9 

2, 3, 6, 9 

2, 3, 6, 9 

7 

7,8 

7,8 

7,8 

7,8,9 

Cs2 

^■1 

- 

- 

- 

HHH 

- 

Cs4 

-13.57 

-12.53 

-11.30 

- r 

-9.92 

-5.21 

C7 

0.09 

0.10 

0.01 

0.01 

0.01 

Cs8 

-0.62 

-0.61 

-0.60 

-0.59 

-0.49 

Cs9 

- 

- 

“ 

-0.01 

0.38 

AF/R 

-6.14  . 

-1.72 

-1.73 

-1.72 

-1.56 

a 

e  -  5  for  q  =  12  -  22 
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TABLE  9.6  (Continued)  DELTA  WING  OPTIMIZATION  WITH 
STRESS  AND  THICKNESS  CONSTRAINTS 


U i 

0.1 

0.1 

0.1 

0.1 

0.1 

M 1 

m 

0.1 

0.1 

_ 

0.1 

0.098 

Design 
SIo.  (q) 

18 

19 

20 

21 

22 

m 

0.2330 

0.1594 

0.1022 

0.1543 

0.1000 

p2 

0.1 

0.1 

0.1 

0.1 

0.1 

P3 

0.1 

0.1 

0.1 

0.1 

0.1 

wm 

0.6248 

0.6120 

0.6603 

0.6123 

0.6589 

m i 

0.1334 

0.1132 

0.1207 

0.1115 

0.1195 

El 

0.1 

0.1 

0.1 

0.1 

0.1 

EB 

1.7754 

1.8032 

1.7417 

1.8039 

1.7440 

p8 

0.9126 

0.9424 

0.9253 

0.9449 

0.9274 

p9 

0.1 

0.1 

0.1 

0.1 

0.1 

PlO 

0.2534 

0.3019 

0.2856 

0.3110 

0.2887 

F(p) 

8.874 

8.807 

8.811 

8.810 

8.80  8 

The. 

Con. 

2, 3, 6, 9 

2, 3, 6, 9 
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F(622)  =  8.808  correspond  to  67.5%  structural  and  54%  total 
weight  reduction. 

In  figure  9.4  the  merit  function  F(p)  is  plotted 
against  step  sizes.  Again  it  can  be  seen  that  the  improvement 
in  the  merit  function  after  11th  design  step  is  poor  compared 
to  the  number  of  design  cycles.  F(PX1)  =  9.884  approximately 
correspond  to  94%  of  the  weight  reduction  achieved  at  22nd  de¬ 
sign  step.  Using  smaller  step  sizes  and  greater  number  of  de¬ 
sign  cycles  best  possible  was  to  obtain  F(p)  =  8.765  with  simi¬ 
lar  thickness  distribution. 

Skin  thickness  contours  corresponding  to  the  minimum- 
weight  design  (£22^  satisfying  stress  and  thickness  constraints 
are  presented  in  figure  9.5. 

9.4.3  Optimization  with  Flutter,  Stress  and  Thickness  Con¬ 
straints 

Finally  the  delta  wing  was  optimized  to  satisfy  the 
flutter,  stress  and  thickness  constraints.  The  flutter  con¬ 
straint  was  considered  active  at  all  design  steps,  stress  and 
thickness  constraints  became  active  as  they  were  violated. 

Table  9.7  presents  the  results  for  18  design  steps. 
These  results  were  obtained  in  three  different  computer  runs, 
six  design  steps  for  each  run.  On  the  CDC  6600  each  design 
cycle  took  approximately  24.5  seconds  which  included  the  calcu¬ 
lation-  of  the  flutter  constraint,  the  flutter  constraint  grad¬ 
ients,  the  critical  stress  constraints,  the  gradients  of  the 
critical  stress  constraints,  critical  thickness  constraints  and 
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TABLE  9.7  DELTA  WING  OPTIMIZATION  WITH  FLUTTER 
STRESS  AND  THICKNESS  CONSTRAINTS 
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TABLE  9.7  (Continued)  DELTA  WING  OPTIMIZATION 
WITH  FLUTTER  STRESS  AND  THICKNESS  CONSTRAINTS 
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TABLE  9.7  (Continued)  DELTA  WING  OPTIMIZATION 
WITH  FLUTTER  STRESS  AND  THICKNESS  CONSTRAINTS 
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their  gradients  and  the  calculation  of  Ap  using  the  gradient 

~q 

projection  algorithm.  In  the  determination  of  active  stress 
constraints,  e  =  25  was  used  for  the  first  run,  e  =  15  was  used 
for  the  second  run  and  e  =  5  was  used  for  the  last  run.  The 
reductions  in  R^,  R^,  R^  are  for  thickness  constraint  vio¬ 
lations  at  nodal  points  9,  5  and  10  respectively.  The  reduc¬ 
tions  in  Rg,  Rg,  Rlg  are  for  stress  constraint  violations  at 
nodal  points  8,  2  and  2  respectively. 

From  table  9.7  it  can  be  seen  that  until  the  12th 
design  step  the  flutter  constraint  is  closely  followed  within 
the  feasible  region.  At  design  step  13,  15  and  18  the  flutter 
occurs  between  second  and  third  modes.  It  is  easier  to  observe 
this  phenomenon  from  table  9 . 8  where  first  three  eigenvalue 
pairs  during  the  optimization  are  given.  As  has  previously 
been  discussed  for  the  optimization  with  only  flutter  and  thick¬ 
ness  constraints?  the  imaginary  parts  of  the  second  and  third 
eigenvalue  pairs  come  closer  during  the  optimization  until  the 
13th  design  step  where  they  are  close  enough  for  the  flutter  to 
occur. 

At  the  17th  design  step  the  merit  function  is 
£■(£17)  =  11.032,  which  is  minimum  among  feasible  designs.  This 
correspond  to  59.2%  structural  and  47.4%  total  weight  reduction. 
This  time  the  merit  function  at  the  11th  design  step, 

F  (£ll)  =  12.442,  corresponds  to  91%  of  the  weight  reduction 
achieved  at  17th  design  step.  Using  smaller  step  sizes  and 
greater  number  of  design  steps,  best  design  obtained  had  a 
merit  function  F(p)  =  10.990.  The  corresponding  thickness  dis- 
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TABLE  9.8  FIRST  THREE  EIGENVALUE  PAIRS 
DURING  OPTIMIZATION 
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•1.0  ^  Fig.  9.6  Progress  of  Delta  Wing  Optimization  with  Flutter 
Stress  and  Thickness  Constraints 
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tribution  was  very  close  to  the  of  table  9.7. 

In  figure  9.6  the  merit  function  F(p)  and  flutter 
constraint  Cf  are  plotted  against  step  sizes.  The  kink  at  the 
merit  function  plot  for  design  number  6  is  due  to  the  reduction 
of  e  from  25  to  15  which  releases  some  of  the  active  stress 
constraints . 

Skin  thickness  contours  corresponding  to  the  minimum- 
weight  design  (P-^)  satisfying  flutter,  stress  and  thickness 
constraints  are  presented  in  figure  9.7. 

9.5  Discussion  of  Results 

The  same  step  size,  R  ,  pattern  was  used  to  obtain 
the  minimum-weight  designs  of  a  cantilevered  delta  wing  with 
three  different  set  of  constraints.  Each  time  more  than  90%  of 
the  total  weight  reduction  was  obtained  within  11  design  cycles. 
It  was  found  out  that  with  a  little  experience  with  a  particu¬ 
lar  optimization  problem;  the  determination  of  a  suitable  step 
size  pattern  is  not  very  difficult.  Experiences  with  different 
step  size  patterns  also  showed  that  the  number  of  design  cycles 
and  the  final  minimum-weight  design  are  not  very  sensitive  to 
small  changes  in  step  size  patterns.  These  results  are  very 
encouraging  as  to  the  capabilities  of  the  multi-constraint 
gradient  projection  algorithm. 

It  was  also  found  out  that  the  minimum-weight  design 
for  stress  and  thickness  requirements  is  radically  different 
from  the  minimum-weight  design  for  flutter  and  thickness  re¬ 
quirements.  This  substantiates  the  results  reported  by  the 
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Langley  Research  Center  (Ref.  29) . 

Experimenting  with  stress  optimization  showed  that 
using  a  large  value  for  e  initially  and  then  gradually  reduc¬ 
ing  it  is  effective  in  reducing  the  number  of  design  steps  re¬ 
quired  to  obtain  the  minimum-weight  design. 

The  observation  of  the  flutter  constraint  revealed 
that  towards  the  end  of  the  optimization  procedure,  flutter 
involving  modes  other  than  the  initial  two  modes  does  occur. 
However  it  also  became  evident  that  any  improvement  in  merit 
function,  after  such  phenomenon  is  encountered,  would  be 
insignificant  unless,  possibly  multiple  flutter  constraints 
were  employed.  This  possibility  has  not  been  investigated 
within  the  scope  of  this  study. 

It  is  interesting  to  note  that  a  design  was  formed  by 
taking  the  larger  of  the  two  nodal  point  thicknesses  between 
the  minimum-weight  design  which  only  satisfy  stress  and  thick¬ 
ness  constraints  and  the  minimum-weight  design  which  only 
satisfy  flutter  and  thickness  constraints.  This  design  had  a 
merit  function  F(p)  =  11.292  and  satisfied  all  stress  con¬ 
straints  but  failed  to  satisfy  the  flutter  constraint. 


CHAPTER  10 


CONCLUSIONS  AND  SUGGESTIONS 

10.1  Conclusions 

The  conclusions  obtained  from  this  study  can  be 
summarized  as  follows: 

1.  The  proposed  procedures  for  handling  the  flutter 
constraint  has  been  proved  to  be  very  satisfactory. 

2.  The  analytical  expressions  for  the  calculation  of 
flutter  and  stress  constraint  gradients  proved  to  be  both 
accurate  and  time  saving. 

3.  The  proposed  optimization  algorithm  which  is  based 
on  gradient  projection  concept,  was  proved  to  be  very  effective 
in  obtaining  approximate  feasible  minimum-weight  designs  within 
few  (10-11)  design  cycles.  It  was  also  found  out  that  the 
efficiency  decreases  as  the  optimum  design  is  approached. 

4.  Experience  with  panel  and  delta  wing  optimization 
problems  suggested  that  establishing  a  step  size  pattern  for  a 
particular  structural  optimization  problem  is  not  a  very  diffi¬ 
cult  task  and  that  once  this  pattern  is  set,  it  can  be  used  to 
optimize  the  same  structural  model  with  different  combinations 
of  constraints.  However,  more  research  is  required  in  order  to 
establish  more  concrete  conclusions. 

5.  Results  of  panel  optimization  indicated  that  the 
effect  of  damping  on  mi numum- weight  designs  was  negligible  for 
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practical  damping  values.  However  results  also  suggested  that 
large  savings  are  possible  for  members  where  a  high  structural 
damping  is  introduced  for  special  reasons. 

6.  Results  of  delta  wing  optimization  showed  that 
optimum  designs  for  stress  and  for  flutter  requirements  are 
radically  different.  This  fact  provides  the  motivation  for 
combined  (aeroelastic-stress)  optimization  approach. 

10 . 2  Suggestions  for  Further  Research 

Several  suggestions  for  the  extension  of  this  work 
are  listed  in  the  following: 

1.  Introduction  of  the  capability  to  handle  multiple 
flutter  constraints:  The  results  of  cantilevered  delta  wing 
optimization  showed  that;  a)  flutter  could  occur  between  modes 
other  than  the  initial  two  modes,  b)  multiple  mode  pairs  could 
be  involved  at  a  particular  design  step.  Such  a  capability 
therefore  could  be  very  effective  in  handling  such  situations 
and  possibly  will  produce  better  designs. 

2.  Introduction  of  a  procedure  to  determine  the  step 
sizes  automatically:  It  is  conceivable  that  a  procedure  can  be 
developed  which  automatically  determines  step  sizes,  possibly 
using  constraint  gradient  information  of  earlier  design  steps. 
However  for  the  present  it  seems  that  the  best  method  to  deter¬ 
mine  the  step  size  pattern  for  an  unfamiliar  structural  opti¬ 
mization  problem  is  to  run  the  optimization  interactively. 

3.  Introduction  of  modal  representation  to  reduce  the 


size  of  the  eigenvalue  problem  which  represents  the  flutter  con¬ 
straint.  Such  a  development  would  involve  the  calculation  of 
required  number  of  natural  modes  and  the  reduction  of  system 
matrices  using  these  modes.  Relative  effectiveness  of,  keeping 
the  initial  mode  shapes  throughout  the  entire  optimization  or 
updating  the  mode  shapes  as  the  design  changes,  should  be 
investigated. 

4.  A  study  should  be  conducted  to  determine  the 
relative  efficiency  of  the  different  optimization  algorithms  on 
aeroelastic-stress  optimization  problems. 


APPENDIX  A 


CONSTANT  THICKNESS  SANDWICH  BEAM  ELEMENT 


Pig.  A.l  Constant  Thickness  Sandwich  Beam  Element 


Figure  A.l  shows  the  cross  section  and  the  degrees  of 
freedom  of  typical  element.  A  design  variable  is  associated 
with  each  element.  We  have 

p(x)  =  p±  (A.l) 

and 

Fe(p)  =  P±  (A. 2) 

where  Fe(p)  is  the  merit  function  associated  with  an  element. 

In  the  following,  non-dimensional  element  stiffness, 
mass,  aerodynamic  and  damping  matrices  are  given  in  the  form 
that  they  were  employed  in  equations  (5.3)  and  (5.4). 
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A.  1  Stiffness  Matrix 
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.2  Mass  Matrix 
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where  n  is  same  as  defined  in  equation  (3.8) 
A. 3  Aerodynamic  Matrix 
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A. 4  Damping  Matrix 
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In  all  the  above  expressions  is 


n 


£ 

a 


(A. 6) 


(A. 7) 


where  £  has  the  same  meaning  as  was  defined  in  Chapter  5,  sec¬ 
tion  5.1.  For  the  panel  problem,  £  was  taken  to  be  the  span 
length. 


APPENDIX  B 


TAPERED  SANDWICH  BEAM  ELEMENT 


Figure  B.l  shows  the  main  properties  of  a  typical 
element.  A  design  variable  pi  is  associated  with  each  nodal 
point  and  we  have 


and 


P  (x) 

"  Pi  +  <pi+l  "  pi>x/a 

(B.l) 

Fe(p) 

■  5(pi  +  pi+l> 

(B .  2) 

Non-dimensional  element  stiffness,  mass,  aerodynamic 
and  damping  matrices  as  they  were  used  in  equations  (5.3)  and 
(5.4)  are  given  in  the  following. 
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B.l  Stiffness  Matrix 


where  n  is  same  as  defined  in  equation  (3.8) 
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B.3  Aerodynamic  Matrix 


(B.  5) 


where  £  is  a  convenient  scaling  dimension  as  was  defined  in 
section  5.1.  For  the  panel  problem  the  span  length  £  was  used 
for  this  purpose. 


APPENDIX  C 


HIGH  PRECISION  TRIANGULAR  SANDWICH  PLATE  BENDING 
ELEMENT  WITH  LINEAR  THICKNESS  VARIATION 

Figure  C.l  shows  the  plan  view  of  a  typical  element. 

It  has  a  core  of  constant  thickness  d  and  cover  skins  with 

linear  thickness  variation.  Each  nodal  point  has  6  degrees  of 

freedom  (W,  W  ,  W  ,  W  ,  VI  ,  W  )  .  The  deflection  W(£,n) 

x  y  xx  xy  yy 

within  a  triangular  element  is  taken  as  a  quintic  polynomial 
(Ref.  22). 

W(5,n)  =  a1  +  a2£  +  a3n  +  a4£2  +  a5£n  +  agr)2  +  a7£3 

t  ag£  h  t  ag£r|  +  a^QTi  ^  t  a^2?  n 

+  al352n2  +  a14?n3  +  a15n4  +  alg£5  +  a17sV 

+  a18^2n3  +  aig$n4  +  a20n5  (C.l) 
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where  R  transforms  nodal  point  degree  of  fredoms  from  local 
coordinates  (£,n)  to  global  coordinates  (x,y),  and  T  transforms 
polynomial  coefficients  to  nodal  point  degree  of  freedoms  in 
local  coordinates.  Details  of  the  derivation  of  T(afbfc)  and 
R(0)  can  be  found  in  reference  22  where  T_1  is  given  in  table 
1  and  R  is  given  in  table  2 . 


where 


Using  equation  (C.5)  for  a  we  can  write, 


W(£,n)  =  9  T  R  We  (C.6) 


tT  =  (K  ln  \  E  2n  2 .  K  20n  20)  (c.7) 


with 

mT  =  {0,1, 0,2, 1,0, 3, 2, 1,0, 4, 3, 2, 1,0, 5, 3, 2, 1,0} 

(C.7a) 

and 

n  =  {0,0, 1,0, 1,2, 0,1, 2, 3, 0,1, 2, 3, 4, 0,2, 3, 4, 5} 

(C.7b) 

which  gives  us  W(5,n)  in  terms  of  nodal  point  degree  of  freedoms 
For  linear  cover  skin  thickness  variation  we  can 
write  p(£,n)  as 


p(CfO)  =  (1  C  n) 


ou 


a. 


a. 


(C .  8) 


From  equation  (C.8)  we  obtain; 


r  \ 


>■  = 


1 

-b 

0  “ 

r  \ 

al 

1 

a 

0 

• 

a2 

1 

0 

c 

a3 

— 

l  J 

(C.  9) 
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or 


*  ■< 

al 

a/r 

b/r 

0 

a2 

= 

-1/r 

1/r 

0 

a3 

*  J 

-a/cr 

-b/cr 

1/c 

r  \ 

pl 

< 

p2  - 

P3 

l  J 

(C.10) 


where  r  =  a  +  b  (Fig.  C.l) .  Thus  we  have 


p  n)  =  U  (  11) 


a/r 

b/r 

0 

r 

pi 

-1/r 

1/r 

0 

« 

P2  ’ 

-a/cr 

-b/cr 

1/c 

p3 

l  °  J 

(C.ll) 


which  gives  us  p(£,ri)  in  terms  of  nodal  point  design  parameters, 
In  the  following  the  element  stiffness,  mass,  aerody¬ 
namic  and  damping  matrices  are  developed  based  on  the  above  re¬ 
lations  and  the  formulas  given  in  reference  22  for  a  constant 
thickness  plate. 


C.l  Stiffness  Matrix 


Using  classical  bending  theory  of  plates,  we  can 
write  the  strain  energy  for  an  element  as 


u 


¥o 


P  (S/n)  {w^  +  w 


■nn  +  2vw«wnn  +  2  Ci-v>w|„}aean 


(C.12) 


where  Dq  is  the  initial  flexural  rigidity  as  discussed  in  Chap¬ 
ter  3,  v  is  the  Poisson's  ratio  and  the  integration  is  over  the 
area  of  the  element. 
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Using  expressions  (C.6)  for  W(£,n)  and  (C.ll)  for 
p(CfTi)  and  carrying  the  integrations  over  the  area  we  obtain 


U 

e 


1 

2 


D 

O 


2S 


k 


TRW 

~  ~  ~e 


(C. 13) 


where 


k 


(C. 14) 


12  3 

The  elements  of  k  ,  k  and  k  are  generated  in  the 
computer  according  to  the  following  expressions: 


Define 


m 

aij 

=  rrunij  (mi  -  1)  (m..  -  1) 

(C .  15) 

n 

aij 

=  n±n  j  (ni  -  1)  (n^  -  1) 

(C. 16) 

mn 
a.  . 
13 

=  2(1  -  v)m.m.n.n.  +  v 

1313 

m. n  .  (m.  - 

131 

1)  (n.  - 

1) 

(C. 17) 

+  V 

m.n. (m.  - 
3  13 

1) (n.  - 

1) 

mij 

=  m±  +  mj 

(C. 18) 

n.  . 
13 

=  ni  +  nj 

(C. 19) 
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Here  itk,  m  j  r  n^,  are  the  same  as  the  exponents  of  equations 
(C. 2)  and  (C.7) . 

Also  define 


Fij  -  aIj  F(mij  -  4'nij>  +  F'”ij-  nij  -  4> 

+  a“  F(mij  - 2-  “i,  -  2> 


(C . 20 ) 


ID 


■  a"j  F(mij  -  3-  nij>  +  a?j  F(mii  +  l-  nii  -  4> 


,  mn  „ ,  - 

+  3ij  F (mi j  -  1,  ntj  -  2) 


(C.21) 


Fij  =  aij  F(mij  *  4'  “ij  +  1)  +  F(l"ij'  nij  -  3) 


mn 


+  aH  F  (m»  .  -  2,  n.  .  -  1) 
13  ij  ID 


(C . 22) 


where  F(m,n)  means  (Ref.  22) 


F(m,n) 


Finally  we  have 


k1. 

ID 


k2. 

ID 

k3  . 
ID 


^mnnd?dn  =  cn+1  (am+1  -  (-b)m+1}  , 

(C .  23) 


(C.24) 

(C . 25) 

(C . 26) 


1 

_2 

a 

F3 

r 

r  •  « 

ID 

cr 

r  1  • 

ID 

1 

F2 

F3 

r 

Fij 

cr 

ID 

If?. 

c  13 


where  a,  b,  c  and  r  are  the  element  dimensions  (Fig.  C.l) . 
Thus  the  element  stiffness  matrix  is 


K 

~e 


=  Dq  RT  TT{Pi  k1  +  p2  k2  +  P3  k3}T  R 


(C .  27) 
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or 


Ke(p)  =  RT  TT  {pjk1  +  p,k2  +  p^k3}  T  R 


2~ 


3* 


(C . 28) 


C.2  Mass  Matrix 

For  sinusoidal  time  dependence,  the  kinetic  energy  of 
an  element  can  be  written  as 


Te  -  5“2 


m(E;,n)  W  d£  dr) 


(C. 29) 


where  a)  is  the  frequency  and  m(£,n)  is 


m(C,n)  =  m0  {p(?,n)hm  +  1  -  nm> 


(C .  30 ) 


as  described  in  Chapter  3.  Here  nm  is  the  same  as  n  (Eq.  3.8) 
of  the  text.  Substituting  the  expression  (C.30)  for  m(£,n)  and 
the  expression  (C.6)  for  W(£,n)  in  equation  (C.29)  and  carrying 
out  the  integration  we  obtain; 


T 

e 

1  2 

2^  mo 

where 

“  =  +  P21 

The  elements 

.  1 
of  m  , 

the  computer  according 

to  the 

Define 

pij  " 

F(mij 

F(mij 

,T  _T_ 


*  ~  ~  ~e 


im 


(C .  31) 


(C .  32) 


(C . 33) 


ID' 


(C . 34) 
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F?  . 
13 


F(mij'  nij  +  D 


(C . 35) 


where  m.^ j ,  n^j  snd  F(m,n)  have  the  same  meaning  as  in  equations 
(C.18),  (C.19)  and  (C.23)  for  the  stiffness  matrix.  In  terms 
of  the  above  definitions  we  have 


1 

m.  .  = 

13 

2 

m.  .  = 

13 

3  1  3 

m.  .  =  —  F .  . 

13  c  13 


m? .  =  F^  . 

13  13 


a 

Fi 

1 

F2 

a 

F3 

r 

Fij 

r 

13 

cr 

13 

b 

F1 

+  i 

F2 

b_ 

F3 

r 

13 

r 

Fij 

cr 

13 

(C . 36) 
(C.37) 
(C. 38) 
(C.39) 


where  a,  b,  c  and  r  are  the  element  dimensions  (Fig.  C.l) . 
Thus  the  element  mass  matrix  is 


~e  ~  mo  ~T  ?T  ^nm'-pl1?1  +  p2?2  +  p3?3;*  +  U  “  nm)  mc}  T  R 


(C.40) 


or 


??e(P)  =  ?T  ?T  ^m^l-1  +  p2~2  +  p3?3^  +  U  ”  nm)  mc}  T  R 


(C . 41) 


C.3  Aerodynamic  Matrix 

We  can  use  equations  (4.10)  and  (4.11)  to  derive  the 
aerodynamic  matrix.  Since  the  air  flow  is  on  both  sides  of  the 
delta  wing,  there  will  be  an  additional  factor  of  2.  Noting 
the  similarity  between  equations  (4.8)  and  (C.6)  we  can  write? 
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u  =  2 


2q 

2  Tv  1/2 


PC  -  1) 


~  lx  ?T  d^dn 


(C.42) 


A 


from  which  we  can  define 


a  = 


*  f  d^dn 


(C .  43) 


A 


For  the  angle  0  (Figure  C.l)  between  the  two  coordi¬ 
nate  systems  we  have 


3  JT  3  ,T  q  3  .T  . 

!  =  31  !  cos  6  "  3n  !  sin  9 


(C . 44) 


Using  equation  (C.7)  for  we  obtain 
_m.+m.-l  n.+n. 


aij  = 


r  -m.+m.-l  n.+n.  Q  rm.  +m .  ja.  +n  .-1  . 

tnu  £  i  3  n  l  j  cos  0-n^Ci  j  n  l  j  sin  0}d£dri 


or 


a.  . 
13 


(C. 45) 


m.  F(m. .  -  1,  n. .)  cos  0  -  n.  F (m . . ,  n. .  -  1)  sin  0 

1  Xj  1J  i  XJ  Ij 

(C . 46) 


where  n^j  and  F(m,n)  have  the  same  meaning  as  in  equations 

(C.18),  (C.19)  and  (C.23)  for  the  stiffness  matrix.  From  the 

above  definitions  we  can  write 

-  T  T 

A  =  2  R  T  a  T  R 

0  ~  ~ 


(C.47) 


It  should  be  noted  that  when  deriving  matrices 
e(p)f  Ae  and  pe  the  element  dimensions  are  scaled  by 

as  discussed  in  Chapter  5.  For  the  delta  wing  problem  the 
length  of  the  cantilevered  edge,  L,  (Fig.  9.1)  is  used  for 
this  purpose. 
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C . 5  Boundary  Conditions 

In  this  section  equations  for  applying  the  free  edge 
natural  (force)  boundary  conditions  on  deflection  curvatures 
are  derived. 


Fig.  C.2  Plate  Free  Edge  and  Free  Edge  Moments 


In  figure  C.2  a  free  edge  that  makes  an  angle  a  with 
the  x-axis  is  shown.  For  the  (n,s)  coordinate  system  the  free 
edge  plate-bending  moments  are  and  Mng .  Free  edge  natural 
boundary  conditions  given  by  Kirchhoff  (Ref.  31)  are 


M 


n 

3M 


Qn  "  3s 


=  0 


ns 


=  0 


(C.53) 

(C.54) 


where  Qn  is  the  free  edge  shear  force.  Equation  (C.54)  requires 
the  calculation  of  transverse  displacement  third  derivatives. 
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For  simplicity  it  was  replaced  in  the  present  optimization 
study  by  the  equation 

Mns  =  0  (C.55) 

which  assumes  Qn  =0.  Using  equation  (3.2)  which  relates  plate¬ 
bending  moments  to  deflection  curvatures  we  can  reduce  equations 
(C . 53)  and  (C.55)  to 


W  +  v  W 
nn  ss 


(C . 56) 


(1  -  v)  Wns  =  0  (C.57) 

where  v  is  the  Poisson's  ratio.  Employing  second  degree  tensor 
transformation  (Ref.  22)  we  can  write 


Wss  =  Wxx  cos2a  +  2Wxy  sin  a  cos  a  +  Wyy  sin2a  (C.58) 

2  9 

Wnn  “  Wxx  sin  a  “  2Wxy  sin  a  cos  a  +  wyy  cos  a  (C.59) 

Wns  =-Wxx  sin  a  cos  01  +  WXy  (cos2a  “  sin2a) 


+  wyy  sin  a  cos  a  (C . 60 ) 

Substituting  these  expressions  for  W  .  W  and  W 

ss  nn  ns 

in  equations  (C.56)  and  (C.57)  we  obtain 

2  2  2  9 

Wxx(sin  a  +  v  cos  a)  +  Wyy(cos  a  +  v  sin  a) 

-  2Wxy  (1  -  v)  sin  a  cos  a  =  0  (C.61) 

Wxx  sin  a  cos  a  ”  Wyy  sin  a  cos  a  “  Wxy(cos2a  -  sin2a)  =  0 


(C . 62) 
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which  are  the  constraint  equations  used  on  the  curvatures  of 
free  edge  nodal  points  (i.e.  nodal  points  3,  6,  10,  9  and  8, 
Fig.  9.1)  for  the  delta  wing  problem. 

It  should  be  noted  that  the  authors  of  reference  22 
do  not  approve  such  attempts  to  satisfy  natural  boundary  con¬ 
ditions  only  at  the  nodal  points  (see  Ref.  22,  section  2.6  for 
a  detailed  discussion) . 


C.6  Merit  Function 


We  can  define 


pe<P> 


p (£,n)d£dn 


(C. 63) 


where  Fg(p)  is  the  merit  function  associated  with  an  element. 
Using  equation  (C.ll)  for  p(£,n)  we  obtain 


Fe(p) 


(Px  +  P2  +  P3) 


(C. 64) 


which  can  be  used  to  derive  the  merit  function  for  the  system 
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